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Abstract 



When hadrons scatter at high energies, strong color fields, whose dynamics is described by 
quantum chromodynamics (QCD), are generated at the interaction point. If one represents these 
fields in terms of partons (quarks and gluons) , the average number densities of the latter saturate 
at ultrahigh energies. At that point, nonlinear effects become predominant in the dynamical 
equations. The hadronic states that one gets in this regime of QCD are generically called "color 
glass condensates". 

Our understanding of scattering in QCD has benefited from recent progress in statistical and 
mathematical physics. The evolution of hadronic scattering amplitudes at fixed impact parameter 
in the regime where nonlinear parton saturation effects become sizable was shown to be similar 
to the time evolution of a system of classical particles undergoing reaction-diffusion processes. 
The dynamics of such a system is essentially governed by equations in the universality class of 
the stochastic Fisher-Kolmogorov-Petrovsky-Piscounov equation, which is a stochastic nonlinear 
partial differential equation. Realizations of that kind of equations (that is, "events" in a particle 
physics language) have the form of noisy traveling waves. Universal properties of the latter can 
be taken over to scattering amplitudes in QCD. 

This review provides an introduction to the basic methods of statistical physics useful in QCD, 
and summarizes the correspondence between these two fields and its theoretical and phenomeno- 
logical implications. 



Resume 

Lors de la diffusion de hadrons a haute energie, d'intenses champs de couleur, dont la dy- 
namique est decrite par la chromodynamique quantique (QCD), sont crees au point d'interaction. 
Si on represente ces champs en termes de partons (quarks et gluons), la densite de ces derniers 
sature a tres haute energie. Les effets non-lineaires deviennent alors dominants dans les equations 
dynamiques. Les etats hadroniques que Ton obtient dans ce regime de la QCD sont generiquement 
appeles "condensat de verre de couleur". 

Notre comprehension de la diffusion en QCD a beneficie de progres recents en physique statis- 
tique et en physique mathematique. On a montre que revolution des amplitudes de diffusion 
hadronique a parametre d'impact fixe dans le regime dans lequel les effets non-lineaires de satu- 
ration des densites de partons deviennent importants est semblable a revolution temporelle d'un 
systeme de particules classiques soumis a des processus de type reaction-diffusion. La dynamique 
d'un tel systeme est essentiellement gouvernee par des equations dans la classe d'universalite 
de l'equation de Fisher-Kolmogorov-Petrovsky-Piscounov stochastique, qui est une equation aux 
derivees partielles stochastique et non-lineaire. Les realisations de telles equations (c'est-a-dire les 
cvcncmcnts, dans un langage de physique des particules) ont la forme d'ondes voyageuses bruitees. 
Les proprietes universelles de celles-ci peuvent etre transposees aux amplitudes d 'interactions en 
QCD. 

Ce memoire est une introduction aux methodes de physique statistique utiles en QCD, et 
resume la correspondance entre ces deux domaines ainsi que ses implications theoriques et phe- 
nomenologiques . 
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Chapter 1 



Introduction to high energy 
scattering in QCD 

What is the origin of the mass of ordinary matter? How are the nucleons "glued" together to form 
stable nuclei? How can we understand the "zoo" of the particles (hadrons) which are sensitive to 
the "strong" force? 

The modern theory of strong interactions, quantum chromodynamics (abreviated QCD; For a 
comprehensive textbook, see Ref. (2]), discovered about 40 years ago, seems to have the ability 
to help all these problems and many others in a most compact and elegant way. This theory is 
parallel in its formulation to the more well-known theory of electromagnetic interactions, quantum 
electrodynamics, and, with some caveats still to be understood, to the theory of weak interactions, 
consistently with the idea of a (partial) unification of the elementary forces. However, QCD poses 
outstanding mathematical problems, and it became soon clear that its various regimes had to 
be explored by dedicated experiments and specialized theoretical tools. While "simple" fixed- 
order perturbation theory has proved extremely successful to investigate electrodynamics due to 
the intrinsic weakness of the force acting between charged leptons (the characteristic coupling is 
a om ~ ) i chromodynamics has to deal with strong coupling instead (the coupling a s is of order 
0.1 in the most favorable cases and up to 1 in general), and with more subtle nonlinearities which 
severely limits the use of fixed-order perturbation theory. Perturbative expansions have to be 
handled with care, and dedicated tools have to be invented for QCD. 

It has happened that methods were borrowed from other fields of physics: For example the 
computation of low-energy properties of the hadrons (masses, decay rates...) are investigated using 
lattice field theory, like in solid state physics. More recently, tools developed by string theorists 
have proved useful to address the calculation of specific processes involving very high-order terms 
in a perturbative expansion. 

This review, which is a revised version of the author's publication (lj, summarizes some recent 
investigations in a specific regime of quantum chromodynamics, the so-called high-energy or, more 
technically, "small- x" regime. We focus on how this regime is formally related to some models which 
appear in statistical physics and show how this correspondence may be used. We shall first go 
over the recent history of high-energy QCD in order to better expose the context of this research. 

Short history of the field. The study of quantum chromodynamics in the high-energy regime 
has undergone a rapid development in the last 15 years with the wealth of experimental data that 
have been collected, first at the electron-proton collider DESY-HERA, and then at the heavy-ion 
collider RHIC. More energy in the collision enables the production of objects of higher mass in the 
final state, and thus the discovery of new particles. But higher energies make it also possible to 
observe more quantum fluctuations of the incoming objects, that is to say, to study more deeply 
the structure of the vacuum. 

Analytical approaches to QCD in this regime are based on a sophisticated handling of per- 
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turbative expansions of observables in powers of the strong coupling constant a s which, thanks 
to asymptotic freedom, is justified for carefully chosen observables in special kinematical regimes. 
Some sophistication is needed because in the evaluation of Feynman graphs, the coupling constant 
always comes with "infrared" and "collinear" logarithms that are related to the phase space that is 
available to the reaction, that is to say, to kinematics, and may easily push the effective coupling 
to large values. Rcsumming part of these logarithms is mandatory. Rcsumming all of them is too 
difficult. The question is to carefully select the dominant ones, and this is not at all easy. 

An experimental facility able to investigate the high-energy regime of QCD was the HERA 
collider, where electrons or positrons scattered off protons at the center-of-mass energy y/s, ex- 
changing a photon of virtuality Q. Through the scattering, one could probe partonic fluctuations 
of the proton (made of quarks and gluons) of transverse momenta k ~ Q, and longitudinal mo- 
mentum fractions x ~ Q 2 /{Q 2 + s). 

For a long time, the dominant paradigm had been that the collinear logarithms InQ 2 , that 
become large when Q 2 is large compared to the QCD confinement scale A 2 , were the most impor- 
tant ones. As a matter of fact, searches for new particles or for exotic physics require to scrutinize 
matter at very small distances, and hence very large Q 2 have to be considered. Perturbative series 
of powers of a s InQ 2 have to be fully resummed. The equation that performs this resummation is 
the celebrated Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation [3j|5]. 

However, once HERA had revealed its ability to get extremely good statistics in a regime in 
which Q 2 is moderate (from 1 to 100 GeV 2 ) and x very small (down to 10 -5 ) it became clear 
that infrared logarithms (lnl/x) could show up and even dominate the measured observables. 
The resummation of the series of infrared logs is performed by the Balitsky-Fadin-Kuraev-Lipatov 
(BFKL) equation [6^]. The series ^(a s In l/x) fe (with appropriate coefficients) is the leading 
order (LO), while the series ^ a s (a s In is the next-to-leading order (NLO), which has also 
been computed |9 10 . The BFKL equation is a linear integro-differential equation. 

At ultrahigh energy, the bare BFKL equation seems to violate the Froissart bound, that states 
that total hadronic cross-sections cannot rise faster than (In 2 s)/m 2 . The latter is a consequence 
of the unitarity of the probability of scattering. The BFKL equation predicts a power rise with the 
energy of the form s £ , where e is positive and quite large (0.3 to 0.5 according to the effective value 
of a s that is chosen). The point at which the BFKL equation breaks down depends on the value 
of the typical transverse momentum which characterizes the observable (It is the photon virtuality 
Q in the case of deep-inelastic scattering). One may define the energy-dependent saturation scale 
Q s (x) in such a way that the BFKL equation holds for Q > Q s (x). For Q ~ Q s (x), the probability 
for scattering to take place is of order 1, and for Q < Q s (x), it would be larger than 1 if one trusted 
the BFKL equation. The saturation scale is a central observable, which we shall keep discussing 
in this review: It signs the point at which the linear (BFKL) formalism has to be corrected for 
nonlinear effects. The regime in which nonlinearities manifest themselves is a regime of strong 
color fields, sometimes called the color glass condensate (For the etymology of this term, see e.g. 



the lectures of Ref. 11 ; for a review, see Ref. 12 ) 



The fact that unitarity is violated is not only due to the lack of a hadronic scale in the BFKL 
equation, which is a perturbative equation; Introducing confinement in the form of a cutoff would 
not help this particular problem: The violation of unitarity which we are talking about occur at 
small distances. It is just that still higher orders are needed. The NLO corrections to the BFKL 
kernel indeed correct this behavior in such a way that the description of the HERA data in the 
small- a; regime is possible by the BFKL equation. However, these corrections are not enough to 
tame the power-like growth of cross-sections as predicted by the LO BFKL equation. It seems 
that a resummation of contributions of arbitrary order would be needed. 

New equations were proposed well before the advent of colliders able to reach this regime. 
Gribov-Levin-Ryskin wrote down a model for the evolution of the hadronic scattering cross-sections 



in the early 80's 1 13 14 , and Mueller and Qiu derived a similar equation from QCD a bit later 
15 . These equations are integral evolution equations with a nonlinear term, which basically 
takes into account parton saturation effects, that is to say, recombination or rescattering. The 
latter cannot be described in a linear framework such as the BFKL formalism. Subsequently, 
more involved QCD evolution equations were derived from different points of view. In the 90's, 
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McLerran and Venugopala n [16 - 18 proposed a first model, mainly designed to approach heavy-ion 
collisions. Later. Balitsky |19 , Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov and Kovner 
(B-JIMWLK) |20}|24] worked out QCD corrections to this model, and got equations that reduce to 
the BFKL equation in the appropriate limit. Technically, these equations actually have the form 
of an infinite hierarchy of coupled integro-differential equations (in Balitsky 's formulation |19| ), of 
a functional renormalization group equation, or alternatively, of a Langevin equation (in Weigert's 



formulation [24]). A much simpler equation was derived in 1996 by Balitsky 19 and rederived 
by Kovchegov in 1999 [25] [26] in a very elegant way within a different formalism. The obtained 
equation is called the Balitsky- Kovchegov equation (BK). The latter derivation was based on 
Mueller's color dipole model 127], which proves particularly suited to represent QCD in the high- 
energy limit. 

The exciting feature of this kinematical regime of hadronic interaction from a theoretical point 
of view is that the color fields are strong, although, at sufficiently high energies, the QCD coupling 
is weak, authorizing a perturbative approach, and thus some of the analytical calculations out- 
lined above. In such strong field regime, nonlinear effects become crucial. But the conditions of 
applicability of the different equations that had been found had never been quite clear. Anyway, 
these equations are extremely difficult to solve, which had probably been the main obstacle to 
more rapid theoretical developments in the field until recently. 

Furthermore, for a long time, the phenomenological need for such a sophisticated formalism was 
not obvious, since linear evolution equations such as the DGLAP equation were able to account for 
almost all available data. But Golec-Biernat and Wiisthoff showed that unitarization effects may 
have already been seen at HERA 1 28 , 29 . Their model predicted, in particular, that the virtual 



photon-proton cross-section should only depend on one single variable t, made of a combination of 
the transverse momentum scale (fixed by the virtuality of the photon Q) and x. This phenomenon 
was called "geometric scaling" |30]. It was found in the HERA data (see Fig. |1.1| : This is maybe 
one of the most spectacular experimental result from HERA in the small- x regime. 

This observation has triggered many phenomenological and theoretical works. Soon after 
its discovery in the data, geometric scaling was shown to be a feature of some solutions of the 
Balitsky-Kovchegov (BK) equation, essentially numerically, with some analytical arguments (see 
e.g. |32j - [35| ) . The energy dependence of the saturation scale was eventually precisely computed by 
Mueller and Triantafyllopoulos [36]. Later, it was shown that the BK equation is actually in the 
universality class of the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation (37 38 , and 
geometric scaling was found to be implied by the fact that the latter equation admits traveling-wave 
solutions [39] . 

A first step beyond the BK equation, in the direction of a full solution to high energy QCD, 
was taken by Mueller and Shoshi in 2004 [40] . Actually, they did not solve the B-JIMWLK 
equations, which would be the natural candidate for a complete theory. Instead, they solved 
the linear BFKL equation with two absorptive boundary conditions, which they argued to be 
appropriate to represent the expected nonlinearities. Geometric scaling violations were found 
from their calculation, which should show up at any energies. 

Subsequently, it was shown that high-energy QCD at fixed coupling is actually in the univer- 
sality class of reaction- diffusion processes, studied in statistical physics, whose dynamics may be 



encoded in equations similar to the stochastic FKPP equation 41 . The Mueller-Shoshi solution 



was shown to be consistent with solutions to the latter equation. So high-energy QCD seems to be 
in correspondence with disordered systems studied in statistical physics. This correspondence has 
provided a new understanding of QCD in the high-energy regime, and it has proven very useful 
to find more features of high-energy scattering. 



Outline of this memoir. The next chapter is devoted to describing scattering in QCD from 
a s-channel point of view, relying essentially on the parton model or, rather, on a realization 
useful in the high-energy limit, the color dipole model. Once this picture is introduced, it is 
not difficult to understand the correspondence with reaction-diffusion processes occuring in one 
spatial dimension, whose dynamics is captured by equations in the universality class of the Fisher- 
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Figure 1.1: [From Ref. (3l|] Photon-proton total cross-section from the most recent set of deep- 
inelastic scattering data in the low-x regime plotted as a function of a single scaling-variable 
t = Q 2 /Q 2 (x), where Q is the virtuality of the photon and Q 2 (x) ~ A 2 x~°' 3 is the so-called 
saturation scale. Although the cross-section is a priori a function of two variables, all data fall on 
the same curve. This phenomenon is called geometric scaling and was discovered in Ref. 1 30 . 
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Kolmogorov-Petrovsky-Piscounov (FKPP) equation. We then explain how traveling waves appear 
in this context. In Chap. [3] we study in greater detail a toy model for which many technics (field 
theory, statistical methods) may be worked out completely. This model however ignores spatial 
dimensions, and thus, does not account for traveling waves. We summarize the state-of-the-art 
research on equations in the universality class of the FKPP equation in Chap. [4] We then come 
back to QCD, discussing the issue of the impact parameter. This will lead us to introduce new 
models beyond the simple one-dimensional reaction-diffusion type models (Chapter [5]) . Finally, 
we will show how noisy traveling waves may show up in the actual data. 

Over the last few years, several hundreds of papers have appeared related to this subject, 
mainly issued from a very active though restricted community. Obviously, this memoir cannot 
give a complete account of this abundant literature. As a matter of fact, some important recent 
developments had to be left out. Concerning the correspondence itself, we do not attempt to 
establish a definite stochastic nonlinear evolution equation for QCD amplitudes, for to our judge- 
ment, this research line is not mature enough yet: A better understanding of the very saturation 
mechanism at work in QCD is definitely needed before one may come to this issue. Furthermore, 
it is not clear to us that a stochastic formulation would be a technical progress, since there are 
not many known methods to handle complicated stochastic equations. We feel that the same is 
true for the search for effective actions that would include so-called Pomeron loops. We also do 
not address the developments based on the boost-invariance symmetry that scattering amplitudes 
should have: This would drive us too far off the main focus of this review. As for more phenomeno- 
logical aspects, we only discuss the basic features of total cross-sections without attempting to 
address other observables such as diffraction. We do also not address the issue of next-to-leading 
effects such as the running of the QCD coupling. This discussion, though crucial if one wants 
to make predictions for actual colliders, would probably only be technical in its nature: There is 
no conceptual difference between the fixed coupling and the running coupling cases. Here, only 
basic phenomenological facts brought about by this new understanding of high-energy QCD are 
addressed, namely geometric scaling and diffusive scaling. 
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Chapter 2 



Hadronic interactions and 
reaction-diffusion processes 

We shall introduce here the physical picture of high-energy scattering in the parton model. In 
the first section, the color dipole model [21^ is described since it is particularly suited to address 
high-energy scattering, especially close to the regime in which nonlinear effects are expected to play 
a significant role. In a second section, we shall argue that high-energy scattering is a peculiar 
reaction- diffusion process. 



Contents 

2.1 Parton model and dipoles 1101 

2.1.1 General picture [TU] 

2.1.2 BFKL equation from the dipole model [11] 

2.1.3 Unitarity and the Balitsky-Kovchegov equation [T3] 

2.1.4 The B-JIMWLK formalism [TBI 

2.1.5 Saturation [17] 

2.1.6 The Pomeron language [18] 

2.2 Analogy with reaction-diffusion processes 1211 

2.2.1 The BK equation and the FKPP equation ED 

2.2.2 Reaction-diffusion processes: an example [2H 

2.2.3 Universality class of high-energy QCD [21] 



2.1 Parton model and dipoles 
2.1.1 General picture 

For defmiteness, let us consider the scattering of a hadronic probe off some given target, in the 
restframe of the probe and at a fixed impact parameter, that is to say, at a fixed distance between 
the probe and the center of the target in the two-dimensional plane transverse to the collision 
axis. In the parton model, the target interacts through one of its quantum fluctuations, made 



of a high-occupancy Fock state if the energy of the reaction is sufficiently high (see Fig. 2.1 1). 
As will be understood below, the probe effectively "counts" the partons in the Fock state of the 
target whose transverse momenta k (or sizes r ~ 1/fc) are of the order of the momentum that 
characterizes the probe: The amplitude for the scattering off this particular partonic configuration 
is proportional to the number of such partons. 

The observable that is maybe the most sensitive to quantum fluctuations of a hadron is the 
cross-section for the interaction of a virtual photon with a hadronic target such as a proton or a 
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nucleus. The virtual photon is emitted by an electron (or a positron). What is interesting with this 
process, called "deep-inelastic scattering", is that the kinematics of the photon is fully controlled by 
the measurement of the scattered electron. The photon can be considered a hadronic object since 
it interacts through its fluctuations into a quark-antiquark state. The latter form a color dipole 
since although both the quark and the antiquark carry color charge, the overall object is color 
neutral due to the color neutrality of the photon. The probability distribution of these fluctuations 
may be computed in quantum electrodynamics (QED). Subsequently, the dipole interacts with the 
target by exchanging gluons. The dipole-target cross-section factorizes at high energy. One typical 
event is depicted in Fig. |2.1| i. 

Dipole models |42[|43| have become more and more popular among phenomenologists since 
knowing the dipole cross-section enables one to compute different kinds of observables. Like 
parton densities, the latter is a universal quantity, that may be extracted from one process and 
used to predict other observables. Different phenomenological models may be tried for the dipole 
cross-section. QCD evolution equations may even be derived, as we shall explain below. A critical 
recent study of the foundations of dipole models may be found in Ref. |44||45| . 

In QCD, the state of a hadronic object, encoded in a set of wave functions, is built up from 
successive splittings of partons starting from the valence structure. This is visible in the example of 
Fig. |2.1^ l: The quark and the antiquark that build up, in this example, the target in its asymptotic 
state each emit a gluon, which themselves emit, later on in the evolution, other gluons. As one 
increases the rapidity y by boosting the target, the opening of the phase space for parton splittings 
makes the probability for high occupation numbers larger. Indeed, the probability to find a gluon 
that carries a fraction z (up to dz) of the momentum of its parent parton (which may be a quark 
or a gluon) is of order a s N c dz/z for small z (N c is the number of colors; N c = 3 in real-life QCD). 
As we see, there is a logarithmic singularity in z. meaning that emissions of very soft gluons (small 
z) are favored if they are allowed by the kinematics. The splitting probability is of order 1 when 
the total rapidity of the scattering y — In 1/x is increased by roughly 1 /a (the convenient notation 
a = a s N c /n has been introduced). Only splittings of a quark or of a gluon into a gluon exhibit the 
l/z singularity. Therefore, at large rapidities, gluons eventually dominate the partonic content of 
the hadrons. 

The parton model in its basic form, where the fundamental objects of the theory (quarks and 
gluons) are directly considered, is not so easy to handle in the high-energy regime. One may 
considerably simplify the problem by going to the limit of a large number of colors (N c 3> 1), in 
which a gluon may be seen as a zero-size quark-antiquark pair. Then, color-neutral objects become 
collections of color dipoles, whose endpoints consist in "half gluons" (see Fig. |2.1| j). There is only 
one type of objects in the theory, dipoles, which simplifies very much the picture. Furthermore, 
going to transverse coordinate space (instead of momentum space, usually used in the DGLAP 
formalism) by trading the transverse momenta of the gluons for the sizes of the dipoles (through an 
appropriate Fourier transform) brings another considerable simplification. Indeed, the splittings 
that contribute to the amplitudes in the high-energy limit are the soft ones, for which the emitted 
gluons take only a small fraction of the momentum of their parent, the latter being very large. 
Therefore, the positions of the gluons, and thus of the edges of the dipoles, in the plane transverse 
to the collision axis are not modified by subsequent evolution once the gluons have been created. 
Thus, the evolution of each dipole proceeds through completely independent splittings to new 
dipoles. 

We will now see how this picture translates into a QCD evolution equation for scattering 
amplitudes, first in the regime in which there are no nonlinear effects. In a second step, we will 
try and understand how to incorporate the latter. 

2.1.2 BFKL equation from the dipole model 

The building up of the states of each hadron is specified by providing the rate at which a dipole 
whose endpoints have transverse coordinates (xq,xi) splits into two dipoles (xq,X2) and (x2,xi) 
as the result of a gluon emission at position xi when the rapidity of the initial dipole is increased. 
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(b) 



Figure 2.1: (a) The scattering of a virtual photon probe off a particular fluctuation of an evolved 
target made of a quark and an antiquark in its bare state. The photon necessarily goes through 
a quark-antiquark pair at high enough energies, when the target is dominated by dense gluon 
states. (What is represented in this figure is actually the inelastic amplitude, which is a cut of the 
total cross-section or of the forward elastic amplitude), (b) In the dipole model, the probe and 
the target may be represented by sets of color dipoles, and the interaction proceeds through gluon 
exchanges. It is now the elastic amplitude that is represented. The curly vertical lines represent 
2-gluon exchanges between pairs of dipoles. 



It is computed in perturbative QCD and reads |27 



dP 
d(ay) 



(x i -> x 02 ,x 12 ) 



\x - xxf 



d 2 x 2 



\Xo - %2 Fl _ x 2\ 27r 



(2.1) 



Thanks in particular to the large- N c limit, dipole splittings are independent. After some rapidity 
evolution starting from a primordial dipole, one gets a chain of dipoles such as the one depicted 
in Fig. [272} 

The elementary scattering amplitude for one projectile dipole (xq, x\) off a target dipole (zq, z±) 
is independent of the rapidity and reads |27j 



T el ((x ,xi), (z ,Zi)) 



n 2 a 2 s 2 \xo ~ ^i| 2 ki ~ zq\ 2 



\x - z \ 2 \xi - zi\ 



(2.2) 



If the target is an evolved state at rapidity y, then it consists instead in a distribution n(y, {zq, z±)) 
of dipoles. The (forward elastic) scattering amplitude A(y, (xq, x\)) is just given by the convolution 
of n and T cl , namely 



A(y, (x ,xi)) 



d 2 z d 2 z\ 
2?r 2tt 



r ol ((x ,xi), (z ,zi))n(y, (z ,zi)). 



(2.3) 



Let us examine the properties of T el . To this aim, it is useful to decompose the coordinates of 
the dipoles in their size vector r a — xq — x\ (resp. = zq — Z\) and impact parameter b a — X °+ Xl 
(resp. bb — z ° J ! 2 Zl ). In the limit in which the relative impact parameters of the dipoles b = b a — bb 
is very large compared to their sizes, we get the simplified expression 



T e \r a ,r b ,b) - a 2 -^ 

V ^ \r a \,\r b \«\b\ S 6 4 



(2.4) 



and thus the scattering amplitude decays fast as a function of the relative impact parameter. If 
instead the relative impact parameter is small (of the order of the size of the smallest dipole) , we 
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Figure 2.2: Schematic picture of a realization of the dipole evolution after the first two steps of 
the evolution ((a) and (b)), and after some larger rapidity evolution (c). In the first step (a), the 
initial dipole (xq, x{) (denoted by a dashed line) splits to the new dipoles (xq, x%) and (x2, x\) (full 
lines). The points represent the edges of each dipole, that is to say, the position of the gluons. In 
the next step (b), the dipole (x2, x\) itself splits in two new dipoles. The splitting process proceeds 
(c) until the maximum rapidity is reached. Many very small dipoles are produced in the vicinity of 
each of these endpoints, due to the infrared singularity visible in Eq. (2.1 1 (Only a fraction of them 
is represented). The zones 1 and 2 in (c), separated by the transverse distance Ab, would evolve 
quasi-independently after the stage depicted in this figure when saturation effects are included 



(See Sec. 5.1 for the corresponding discussion). 
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get for disymmetric sizes 



T c \r a ,r b ,b) 



r 2 

8 2 ' 



(2.5) 



where r < = min(|r | 
performed. 



= max(|r a |, and where the integration over the angles has been 



Equation (2.4 1 means that the dipole interaction is local in impact parameter: It vanishes as 



soon as the relative distance of the dipoles is a few steps in units of their size. Eq. (2.5 1 shows 



that only dipoles whose sizes are of the same order of magnitude interact. These properties are 



natural in quantum mechanics. Thus the amplitude A in Eq. (|2.3| effectively "counts" the dipoles 

- Xi|), with a weight 



of size of the order of Foil at the impact parameter 
factor a 2 



Xg+xi 



(up to |x i| 



Fo ■ 



An evolution equation for the amplitude A with the rapidity of the scattering can be estab- 
lished. It is enough to know how the dipole density in the target evolves when rapidity is increased, 
since all the rapidity dependence is contained in n in the factorization ( |2.3[), an d such an equation 
may easily be worked out with the help of the splitting rate distribution (2.1 1. It reads \27 



dn(y, (x ,xi)) 
d(ay) 



d 2 x 2 



Foil 



2tt 



|x 2| 2 |xi2| 2 



[n(y, (x , x 2 )) + n(y, (x 2 , ii)) - n(y, (x , xi))], (2.6) 



where x a j, = x a — x b . The very same equation holds for A. The elementary s catt ering amplitude 
T el only appears in the initial condition at y — 0, which is not shown in Eq. (2.6). 



In a nutshell, the integral kernel encodes the branching diffusion of the dipoles. The total 
number of dipoles at a given impact parameter grows exponentially, and their sizes diffuse. The 
appropriate variable in which diffusion takes place is ln(l/|xoi| 2 ). (This is due to the collinear 
singularities in Eq. (2.1).) Equation (2.6 1 is nothing but the BFKL equation. A complete solution 
to this equation, including the impact-parameter dependence, is known |46| . 

An important property of the amplitude A is that it is boost-invariant. This property is 
preserved in the BFKL formulation. We could have put the evolution in the projectile instead 
of the target, or shared it between the projectile and the target: The result for the scattering 
amplitude would have been the same. In a frame in which the target carries y' units of rapidity 
and the projectile y — y', the amplitude A reads 



n 



A(y, (x ,xt)) 



projectile ( q . 



d 2 z d 2 Z\ drz' d 2 z[ 
2ir 2ir 2tt 2tt '' 



projectile 



(y - y', (zo,zi)\(x ,xx)) 

x T c \(z , zi), (4, ^)K argct (y', (4 *!))■ (2-7) 



(y — y', (z , Zi)\(xq,Xx)) is the density of dipoles (z , Zi) found in a dipole of initial size 



(xo,xi) after evolution over y — y' steps in rapidity. If y' = y, one recovers Eq. (2.3 1. If y' = 0, 
then all the evolution is in the projectile instead. 

The amplitude A is related to an interaction probability, and thus, it must be bounded: In 
appropriate normalizations, A has to range between and 1. But as stated above, the BFKL 
equation predicts an exponential rise of A with the rapidity for any dipole size, which at large ra- 
pidities eventually violates unitarity. Hence the BFKL equation is not the ultimate representation 
of high-energy scattering in QCD. 



2.1.3 Unitarity and the Balitsky-Kovchegov equation 

It is clear that one important ingredient that has been left out in the derivation of the BFKL 
equation is the possibility of mult iple scatterings between the probe and the target. Several 
among the n P r °j cctlle dipoles in Eq. (2.7) may actually interact with the n target dipoles in the other 
hadron simultaneously. The only reason why such interactions may not take place is that T cl ~ a 2 
(see Eq. (2.2 1), and thus the probability for two simultaneous scatterings is of order oq , which is 
parametrically suppressed. But this argument holds only as long as the dipole number densities 
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Figure 2.3: Derivation of the Balitsky-Kovchegov equation. 



are of order 1. If n ~ l/a:^ (which is also the point above which the unitarity of A is no longer 
preserved in the BFKL approach), then it is clear that multiple scatterings should occur. 

In order to try and implement these multiple scatterings, we introduce the probability that 
there be no scattering between a dipole (xq, x{) and a given realization of the target in a scattering 
with total rapidity y, that we shall denote by S(y, (xq, Xi)). Let us start with a system in which the 
evolution is fully contained in the target. We increase the total rapidity by boosting the projectile 
(initially at rest) by a small amount dy. Then there are two cases to distinguish, depending on 
whether the dipole (xo,xi) splits in the rapidity interval dy. In case it splits into two dipoles 
(xo, X2) and (x 2 , x±), the probability that the projectile does not interact is just the product of the 
probabilities that each of these new dipoles do not interact. This is because once created, dipoles 
are assumed to be independent. In summary: 



■ S(y, (in, x-i )) with proba 1 — ady f ,f p - s (xm — > xm, x-12) 
S(y + dy,(x ,x 1 )) = { J' 1 1 t -f °p < ^ ^ 

^b{y, {xo,x 2 ))S{y, (x 2 ,xi)) with proba Cfdyj^(x i ->• £02,^12) 



Taking the average over the realizations of the target and the limit dy — ► 0, we get 
d f d^x^ x^ 

■K-{S(y, {x ,x x ))) = a / — j^[(S(y, (x ,x 2 ))S(y, (x 2 ,xi))) - (S(y, (x ,x x )))] (2.9) 

ay j z7r Xq2 x 2 i 



(See Fig. 2.3 for a graphical representation.) We see that this equation is not closed: An evolution 
equation for the correlator (S(y, (xq, x 2 ))S(y, (x 2l Xi))} is required. However, we may assume that 
such correlators factorize in the following sense: 

(S(y,(x ,x 2 ))S(y,(x 2 ,x 1 ))) = (S(y,{x ,x 2 )))(S(y,(x 2 ,x 1 ))). (2.10) 

This assumption is justified if the dipoles scatter off uncorrelated targets, for example, off different 
nucleons of a very large nucleus. Writing A = 1 — (S), we get the following closed equation for A: 

—A(y,(x ,xi)) = a / — 3— ~-\A(y, (x ,x 2 )) + A(y, {x 2 ,x 1 )) 

dy J 2tt xfexit 

- A(y> ( x o,xi)) - A(y, (x ,x 2 ))A(y, (x 2 ,x 1 ))}, (2.11) 
which is the Balitsky-Kovchegov (BK) equation |25|[26]. Note that if one neglects the nonlinear 



term, one gets back the BFKL equation (2.6 1 (written for A instead of n). A graphical represen- 
tation of this equation is given in Fig. |2.4] 

It is not difficult to see analytically that the BK equation preserves the unitarity of A: When 
A becomes of the order of 1, then the nonlinear term gets comparable to the linear terms in 
magnitude, and slows down the evolution of A with y, which otherwise would be exponential. 
Hence the solution of the BK equation will exhibit essentially two regimes: A BFKL regime of 
low density in which A « 1 and in which the evolution proceeds linearly, and a high-density 
regime A ~ 1. At fixed y, the transition between these two regimes is controlled by the so-called 
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Figure 2.4: Picture of the BK equation. All the QCD evolution is put in the probe, which carries 
the total rapidity. It develops a high occupancy state of dipoles, which scatter independently off 
the target. 



S(y,(x ,x 2 )) 

S(yX^2 ,x 3 )) 
S(y,(x, ,Xj )) 



S(y,(x ,x 



» 



X, 



X, 



(a) 



(b) 



Figure 2.5: (a) Contribution to the B-JIMWLK equation for the 2-point correlator restricted to 
dipoles (x' 2 is taken equal to x 2 in this figure), (b) A graph that would also contribute to the 
2-point correlator and that is missing in the B-JIMWLK formalism. 



saturation scale Q s {y), which is the inverse size of the dipoles which scatter with an amplitude A 
equal to some fixed number of order 1, for example \. (Q s also depends a priori on the impact 
parameter, as will be discussed in Chap. [5]). 

Let us go back to Eqs. (2.8 1,(2.9 1 and instead of assuming the factorization of the correla- 
tors (2.10), work out an equation for the two-point correlator (SS). From the same calculation as 
before, we get 

a 



dy 



2'1 



cPxz x 



02 



27T x\^x\ 2 



((S03S32S: 



2>l) 



- (Sq 2 S 2 >1)) 

d?x x^ 

2 12 2 (( , S , 2'3'S'3i<S'o2) — (SwSy-i)) , (2.12) 

a '13 a '32' 



where we have introduced the notation S a b = S(y, (x a ,Xb))- (See Fig. 2.5 1 for the corresponding 
graphical representation.) 

This equation calls for a new equation for the 3-point correlators, and so on. The obtained 
hierarchy is nothing but the Balitsky hierarchy [19] (see also Ref. 47-49]) restricted to dipoles. 



2.1.4 The B-JIMWLK formalism 

For completeness, let us briefly mention the Balitsky-Jalilian Marian-Iancu-McLerran-Weigert- 
Leonidov-Kovner (B-JIMWLK) formalism. It is a systematic approach beyond the large- N c limit, 
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initially supposed to completely describe high-energy scattering under appropriate well-defined 
approximations. (Nowadays, the status of this approach is less clear, as we shall discuss later on). 
Instead of evolving dipole or parton number densities, B-JIMWLK evolves observables constructed 
with the help of Wilson operators 



U(x) = Pexp ig s / cte'M"(x)i a 



(2.13) 



where is the color 4-potential over which one eventually averages to arrive at observables, t a 
the standard color matrics (fundamental representation), and g s = ^JAitas the coupling of the line 
to the color field. The integration in the exponent goes along the trajectory of the quark whose 
interactions are described by U, which is essentially a straight lightlike line at very high energies. 

There are several equivalent ways to present this physics. Let us exhibit the so-called "Balitsky 
hierarchy", which is an infinite system of coupled integro-diffcrential equations for the correlators 
of U. In terms of the C/'s, the S-matrix for dipole scattering off a particular field configuration 
reads 

S ab = ^-TrU{x a )U\x b ). 



(2.14) 



The observables are the associated correlators: For example the dipole amplitude is obtained from 
S by averaging over the color fields. 



The first equation of the Balitsky hierarchy for the observable (Sqi) is identical to Eq. (2.9 1 
The second equation, needed to solve the first one, reads 



dy 



1 

N 2 



21/ 



d 2 x 3 
2?r 



' 02 



2 2 
X 03 X 32 



((/S03S32 — Sq 2 )S 2 i) + 



(go - X 3 ) ■ (Xi - X 3 ) 
(x ~ X 3 ) 2 (X! - X 3 ) 2 



1 



^03^31 
(X - X 3 ) ■ (x 2 ~ X 3 ) 



2 01 2 ('S'02('S'03'5'31 — ^0l)) 

Xr\r> X '. 



(x 2 - x 3 ) ■ (xi - x 3 ) 



(x 2 -.x 3 ) 2 (x - x 3 ) 2 (x 2 - x 3 ) 2 (x 2 - x 3 ) 2 (xi - x 3 ) 2 _ 

(2.15) 



±r (Tr {U X0 UlU X3 UlU X2 Ul)) + ]j- (Tr (U Xo Ul s U X2 Ut x U x M 2 )) - ^01 



The first line is the same as Eq. (2.12). The other terms involve sextupoles and are suppressed at 



large- iV c . Hence in the dipole approximation, we recover Eq. (2.12). 

The factorized correlators (2.101 is a solution of the whole dipole hierarchy, and turn actually 



out to be a good approximation to the solution of the full B-JIMWLK equations. This statement 
was first made after the numerical solution to a version of the B-JIMWLK formalism was worked 
out in Ref. |50| . We note however that the latter simulations did not cover a very large range 
in rapidity, and therefore, they may have missed physical effects that would differentiate the full 
B-JIMWLK equation from its approximate forms. 

We may wonder why there are no terms involving one-point functions in the right-hand side of 
the previous equation. Actually, such terms would correspond to graphs like the one of Fig. |2.5b , 
in which, for example, two dipoles merge. They are expected to occur if saturation is properly 
taken into account. While the restriction of the Balitsky equation to dipoles does a priori not 
drastically change the solution for the scattering amplitudes, such terms would instead have a 
large effect, as we shall discover in the next section. To simplify the discussion, we will stick to 



the dipole approximation, which leads to the evolution equations (2.9 1,(2.12 1 



2.1.5 Saturation 

The BK equation may be well-suited for the ideal case in which the target is a nucleus made of 
an infinity of independent nucleons. But it is not quite relevant to describe the scattering of more 
elementary objects such as two dipoles (or two virtual photons, to be more physical). 

Indeed, following Chen and Mueller (5l] (see also Ref. [52]), let us consider dipole-dipole 
scattering in the center-of-mass frame, where the rapidity evolution is equally shared between the 
projectile and the target (see Fig. 2.6 1). Then at the time of the interaction, the targets are 



17 



CHAPTER 2. HADRONIC INTERACTIONS AND REACTION-DIFFUSION PROCESSES 



dipoles that stem from the branching of a unique primordial dipole. Obviously, the assumption of 



statistical independence of the diffusion centers, which was needed for the factorization (2.101 to 
hold, is no longer justified. 



So far, we have seen that nonlinear effects which go beyond the factorization formula (2.7) are 
necessary to preserve unitarity as soon as n ~ l/ a s- This came out of an analysis of Eq. (2.7 1 
in the restframe of the target. The rapidity j/bfkl at which the system reaches this number of 
dipoles and hence at which the BFKL approach breaks down may be found from the form of the 
typical growth of n with y, namely n{y) ~ e ay . Parametrically, 

2/bfkl ~ 4 In \. (2-16) 



Now we may go to the center-of-mass frame, where Eq. (2.7l with y' = y/2 would describe the 
scattering amplitude in the absence of nonlinear effects. There, the typical number of dipoles in 
the projectile and in the target are well below 1/a 2 : n(j/BFKL/2) ~ l/a s . We actually see that the 
evolution of the dipoles in each of theses systems remains linear until y = 2j/bfkl- In that rapidity 
interval, nonlinear effects consist in the simultaneous scatterings of several dipoles from the target 
and the projectile but the evolution of n still obeys the BFKL equation. Now, performing a boost 



to the projectile restframe, the evolution goes into the target. Formula (2.3 1 should then apply 
for the amplitude A. But if the evolution of the target were kept linear, then the amplitude would 
break unitarity because the number of dipoles would be larger than 1/a 2 . Hence, through some 
nonlinear mechanism, which was represented by multiple scatterings between linearly evolving 
objets in the center-of-mass frame, the dipole number density has to be kept effectively lower 
than I /a 2 , in order to preserve unitarity. This is called partem saturation. The precise saturation 
mechanism has not been formulated in QCD. It could be dipole recombinations due to gluon 



fusion, multiple scatterings inside the target which slow down the production of new dipoles 53 



"dipole swing" as was proposed more recently 54 55 , or any other mechanism. Some of these 



mechanisms may be implemented in simplified toy models; see Chap. [3] 

Hence, unitarity of the scattering amplitudes together with boost-invariance seem to require 
some sort of saturation of the density of partons. It is not clear whether such a mechanism is 
included in the B-JIMWLK formalism, since the latter is not obviously boost-invariant. What is 
clear is that some saturation mechanism lacks in the dipole model. 

A pedagogical review of saturation and the discussion of the relationship between saturation 
and unitarity may be found in Ref. 56 . Original papers include Refs. (57 58 . 



Visualizing saturation: Evolution in different models 

We now wish to illustrate how the different schemes of unitarization (BK unitarization, multiple 
scatterings in the center-of-mass frame, explicit parton number saturation) affect the evolution of 
scattering amplitudes. 

In Fig. |2.7| we plot the 5*-matrix element at different rapidity and as a function of the (loga- 
rithmic) dipole size resulting from the evolution of toy models with dynamics similar to QCD (The 
use of such models will be justified later, when we will establish the correspondence of high-energy 
QCD with more general processes). We see that S goes to zero in a region of sizes that extends 
with rapidity. This phenomenon is slower when there are saturation effects explicitly included in 
the evolution, as discussed above in this section. 



2.1.6 The Pomeron language 

So far, we have presented in detail a s-channel picture of hadronic interactions, and it is in this 
formalism that we will understand most easily the link with reaction-diffusion processes. In the 
s-channel formulation, all the QCD evolution happens in the form of quantum fluctuations of 
the interacting hadrons. However, a picture maybe more familiar to the reader belonging to the 
"traditional" QCD community is a i-channel picture, where the rapidity evolution is put in the 
t-channel, while the projectile and target stay in their bare states. This picture directly stems 
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(a) (b) 

Figure 2.6: (a) Scattering in the dipole model in the center-of-mass frame. The evolution is 
shared between the target and the probe. The amplitude is unitarized through the multiple 
scatterings occurring between the two evolved wave-functions, (b) Boost of the previous graph 
to the restframe of the projectile. There is now twice as much evolution in the target and the 
nonlinear effects should occur inside its wave-function, in the course of the evolution. They may 
take the form of "internal" rescatterings (as depicted), or dipole merging... 




Figure 2.7: [From Ref. [59]] 5-matrix element as a function of the logarithm of the dipole size 
for y = 5,10,20,30 (from the center of the figure towards the outskirts). 1/a 2 = 20. For low 
rapidities (y = 5 and y = 10), the evolution is linear (of BFKL type). At y = 20, the unitarity 
limit has been reached in all calculations (S becomes zero in some region of sizes). Later (y — 30), 
the region in which S — propagates outwards. Different unitarization mechanisms are tested: 
simple BK unitarization ("Unitarized-LAB") , center-of-mass unitarization ("COM") and intrinsic 
saturation. 
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Figure 2.8: The BFKL Pomeron is a sum of ^-channel gluon Feynman diagrams. 




(a) 



(b) 



Figure 2.9: (a) Example of a diagram contributing to the BK equation in the t-channel represen- 
tation (see Fig. 2.4 1. The dashed lines represent Pomerons. The rapidity is proportional to the 
length of the Pomeron lines in the i-channel. (b) Pomeron representation of a class of diagrams 
to which Fig. |2.6^t belongs. 



from the usual Lorentz-invariant formulation of quantum field theory, while the dipole picture (or 
the parton model) is derived in the framework of time-ordered perturbation theory. 

Both pictures have their respective advantages and drawbacks. The covariant formulation 
seems to be more suited for higher-order systematic calculations, since for a given observable the 
number of diagrams is smaller than in the time-ordered (s-channel) formalism. The time-ordered 
formalism seems unpractical beyond the tree-level approximation. On the other hand, the latter 
gives maybe a more intuitive picture of scattering processes and seems to be particularly useful to 
study the approach to unitarity. 

In the i-channel picture, classes of Feynam diagrams can be grouped into "Pomerons" (or 
Reggeized gluons, see Fig. 2.8 1, in terms of which scattering processes may be analyzed. (A 



pedagogical review on how to derive the BFKL equation in such a formalism is available from 
Ref. 1 60 ) . An effective action containing Pomeron fields and Pomeron vertices may be constructed. 
In these terms, the s-channel diagrams of Fig. |2.4| and |2.6^ may be translated in terms of the 
diagrams of Fig. 2.9 The effective action formalism was initially developped in Refs. 61-63 . 
More recently, there has been some progress in the definition of the effective action [64] , some of 



it with the help of the correspondence with statistical physics processes 65 66 
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We will not expand on this formulation in the present review, because it is difficult to see the 
analogy with statistical physics in this framework. A s-channel picture is much more natural. 
However, a full solution of high-energy QCD may require to go back to that kind of calculation 
and compute accurately the 1 — > n Pomeron vertices. This program was formulated some time 
ago 67,68 , and there is continuing progress in this direction (see e.g. [69||70]). 



2.2 Analogy with reaction-diffusion processes 

We are now in position to draw the relationship between high-energy QCD and reaction-diffusion 
processes. In the first section below, we will show that the BK equation is, in some limit, an 
equation that also appears in the context of statistical physics. Second, we will exhibit a simple 
reaction-diffusion model, and show in the final section how this model is related to scattering in 
QCD, even beyond the approximations implied in the BK equation. 



2.2.1 The BK equation and the FKPP equation 

Let us first show at the technical level that under some well-controlled approximations, the BK 
equation (2.11 1 may be mapped exactly to a parabolic nonlinear partial differential equation. This 
observation was first made in Ref. [39] . 

To simplify, we will look for impact-parameter independent solutions: A(y, (xq, x\)) is assumed 
to depend on y and xoi only, not on X °+ Xl , We switch to momentum space through the Fourier 
transformation 

A{y,k)= [ p^-e^ A(y, x i). (2.17) 
This transformation greatly simplifies the BK equation 25 , 26 . It now reads 



d &y A{y, k) = X (-d lnk 2)A(y, k) - A\y, k). (2.18) 

The first term on the right-hand side, which is a linear term, is actually an integral transform 
whose kernel, obtained by F ourier transformation of the BFKL kernel (first three terms on the 
right-hand side of Eq. (2.11)). It is most easily expressed in Mellin space since the powers fc~ 27 
are its eigenfunctions, with the corresponding eigenvalues 

X ( 7 ) = 2^(1) -V>(7)" ^(1-7). (2-19) 

This kernel may be expanded around some real 7 = 70, fixed between and 1. Keeping the terms 
up to C(( 7 — 7 o) 2 ) is the well-known diffusive approximation, which is a good approximation at 
large rapidities. Introducing the notations xo = x(7o), Xo = X'(7o) an d Xo = x"(7o)> the BK 
equation reads 

d &v A = 4dlk^ + (70X0 - Xo)dink*A + (xo - 7oX + 2 ¥ { )^ - A 2 . (2.20) 
Through the linear change of variable (ay,lnfc 2 ) — > (t, x), 

t 



ay = 

Xo - loXo 

Ink 2 = 



Xo _ 7oXo - Xo 



(2.21) 



2(xo - 7oX ) + 7oXo Xo - 7oX + ^ 



one may get rid of the first-order partial derivative in the right handside. We then find that the 
rescaled function 

u{tjX) ^A(y(t)Mk 2 (t,x)) 
Xo - 7oX + 
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obeys the equation 

du(t,x) d 2 u(t,x) 
dt dx 2 



+ u(t,x) - u 2 (t,x), (2.23) 



which is the Fisher [37] and Kolmogorov-Petrovsky-Piscounov [38] (FKPP) equation. This equa- 
tion was first written down as a model for gene propagation in a population in the limit of large 
number of individuals. But it turns out to apply directly or indirectly to many different physical 
situations, such as reaction-diffusion processes, but also directed percolation, and even mean-field 
spin glasses |71| . A recent comprehensive review on the known mathematics and the phenomeno- 
logical implications of the FKPP equation can be found in Ref. [72 . 

As a side remark, we note that if 70 is chosen such that the equation x(lo) — 7ox'(7o) is 
verified, then the mapping drastically simplifies. Actually, this choice has a physical meaning, as 
we will discover in Chap. [4] when we try and solve the BK equation. 



Beyond the exact mapping (2.23) between an approximate form of the BK equation and the 
FKPP equation, the full BK equation is said to be in the universality class of the FKPP equation. 
All equations in this universality class share some common properties, as will be understood below. 
The exact form of the equation is unessential. As a matter of fact, recently, it has been checked 
explicitly that the BFKL equation with next-to-leading order contributions to the linear evolution 
kernel (but keeping the QCD coupling fixed) is also in the same universality class. A mapping to a 
partial differential equation (which involves higher-order derivatives in the rapidity variable) was 
exhibited [73] . What defines physically the universality class of the FKPP equation is a branching 
diffusion process with some saturation mechanism. The details seem unimportant. 

We must however keep in mind that there is for the time being no theorem that would clearly 
state the necessary and sufficient conditions for a model to belong (or not) to the universality class 
of the FKPP equation: Our statements are nothing but conjectures, supported by arguments and 
checked against numerical simulations. 

In the next section, we shall give a concrete example of a reaction-diffusion process: We will 
see how the FKPP equation appears as a fluctuationless (or "mean-field") limit of some stochastic 
reaction-diffusion process. In Ref. |39 , it had not been realized that the analogy of QCD with 
such processes is in fact much deeper than the formal mapping between the BK equation and the 
FKPP equation that we have just outlined. But this is actually the case, as we shall shortly argue. 



2.2.2 Reaction-diffusion processes: an example 

We consider the reaction-diffusion model which was introduced in Ref. |74| . It consists in a set 
of particles which are evolving in discrete time on a one-dimensional lattice. The following rules 
define the dynamics of the system: At each timestep, a particle may jump to the nearest position 
on the left or on the right with respective probabilities pi and p r , and may split into two particles 
with probability A. We also allow each of the n(t,x) particles on site x at time t to die with 
probability Xn(t,x)/N. 

We can guess what a realization of this evolution may look like at large times. The particles 
branch and diffuse (they undergo an evolution which can be represented by a linear finite difference 
equation) until their number n becomes of the order of N, at which point the probability that 
they "die" starts to be sizable, in such a way that their number never exceeds iV by a large 
amount, on any site. If the initial condition is spread on a finite number of lattice sites, the linear 
branching-diffusion process may always proceed towards larger values of |x|, where there were no 
particles in the beginning of the evolution. Hence after some lapse of time (typically larger than 
In TV) a realization will look like a double front connecting an ensemble of lattice sites where a 
quasi-stationary state in which the number of particles is N (up to fluctuations) has been reached, 
to an ensemble of empty sites. One front will move towards x — > +00 , the other one towards 
x — > — 00 as the branching diffusion process proceeds. Let us focus on the front traveling to the 
right. The position of the front X(t) may be defined in different ways, leading asymptotically to 
equivalent determinations, up to a constant. For example, one may define X(t) as the rightmost 
bin in which there are more than N/2 particles, or, alternatively, as the total number of particles 
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n(t,x) n(t+dt,x) 




ooo«ooo« 
oo»ooo«« 




Figure 2.10: Picture of a realization of the system of particles at two successive times. In the 
bins in which the number of particles is of order N, some particles disappear, others are created 
by splittings, but overall the number of particles is conserved up to fluctuations of order \/~N. In 
the bins in which n is small compared to N, the dynamics is driven by branching diffusion. As a 
result, n(t,x) looks like a noisy wave front moving to the right. 



in the realization whose positions are greater than 0, scaled by 1/N. A realization and its time 
evolution is sketched in Fig. |2.10| 

Between times t and t + At, ni(t,x) particles out of n(t,x) move to the left and n r (t,x) of 
them move to the right. Furthermore, n + {t,x) particles are replaced by their two offspring at x, 
and n_(t,x) particles disappear. Hence the total variation in the number of particles on site x 
reads 



n{t + At, x) — n(t, x) = —ni(t, x) — n r (t, x) — n_(i, x) 

+ n+(t,x) + ni(t,x + Ax) + n r (t,x- Ax). (2.24a) 

The numbers describing a timestep at position x have a multinomial distribution: 

n! 



P({m,n r ,n + ,n-}) 



VTVV A™ + (An/TV)"- (l—pi—p r — \— Xn/N) An , (2.24b) 



ni!n r !n + !n_!An! 



where An = n — rii — n r — n + — n_ , and all quantities in the previous equation are understood at 
site x and time t. The evolution of u = n/N is obviously stochastic. One could write the following 
equation: 



i(t + At, x) = (u(t+At, x)) + y/(u 2 (t+At,x)) - (u(t+At,x)) 2 v(t + At, x) 



(2.25) 



where the averages are performed over the time step that takes the system from t to t + At. They 
are conditioned to the value of u at time t. v is a noise, i.e. a random function. The equation was 
written in such a way that it has zero mean and unit variance. Note that the noise is updated at 
time t + At in this equation. 

One can compute the mean evolution of u = n/N in one step of time which appears in the 
right-hand side of Eq. ( |2.25[ ) from Eq. ( |2.24[ ). It reads 

(u(t+At,x)\{u(t,x)})=u(t,x)+pi [u(t,x+Ax)-u(t,x)} 

+p r [u(t,x-Ax)-u(t,x)]+\u(t,x)[l-u(t,x)}. (2.26) 



The mean evolution of the variance of u that appears in Eq. (2.251 may also be computed. The 
precise form of the result is more complicated, but roughly speaking, the variance of u after 
evolution over a unit of time is of the order of u/N for small u ~ 1/N. This is related to the fact 
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that the noise has a statistical origin: Having n particles on the average in a system means that 
each realization typically consists in n ± s/n particles. 



When N is infinitely large, one can replace the u's in Eq. (2.261 by their averages: This would 
be a mean-field approximation. Obviously, the noise term drops out, and the equation becomes 
deterministic. Note that if we appropriately take the limits Ax — > and At — > 0, setting 

A = At, Pr.=Pl = tP^, (2.27) 
(Ax)^ 



the obtained mean-field equation is nothing but the FKPP equation (2.231. For the numerical 
simulations of this model that we will perform in Sec. [4] we will keep At and Ax fixed, which is 
more convenient for computer implementation. 

Thus we have seen that the evolution of reaction-diffusion systems is governed by a stochastic 



equation (2.251 whose continuous limit (At — > 0, Ax — > 0) and mean-field limit (N ^> 1) is a 
partial differential equation of the form of (exactly actually, in our simple case study) the FKPP 
equation. We shall now argue that partons in high-energy QCD form a similar system. 

2.2.3 Universality class of high-energy QCD 

Let us come back to the QCD dipole model discussed in Sec. |2.1| We have seen that evolution 
proceeds through a branching diffusion process of dipoles. Let us denote by T(y, r) the scattering 
amplitude of the probe dipole off one particular realization of the target at rapidity y and at a given 
fixed impact parameter. This means that we imagine for a while that we may freeze the target 
in one particular realization after the rapidity evolution y, and probe the latter with projectiles 
of all possible sizes. Of course, this is not doable in an actual experiment, not even in principle. 
But it is very important for the statistical picture to decompose the physical observables with the 
help of such a "gedanken observable". The amplitude A, which is related to the measurable total 
cross-section, is nothing but the average of T over all possible realizations of the fluctuations of 
the target, namely 

A(y,r) = (T(y,r)). (2.28) 

The branching diffusion of the dipoles essentially occurs in the ln(l/r 2 ) variable. The scattering 
amplitude is roughly equal to the number of dipoles in a given bin of (logarithmic) dipole size, 
multiplied by a 2 s . From unitarity arguments and consistency with boost-invariance, we have seen 
that the branching diffusion process should (at least) slow down in a given bin as soon as the 
number of objects in that very bin is of the order of N = 1/a 2 , in such a way that effectively, the 
number of dipoles in each bin is limited to N. A typical realization of T is sketched in Fig. |2.11| 
As in the case of the reaction-diffusion process, from similar arguments, it necessarily looks like a 
front. The position of the front, defined to be the value r s of r for which T is equal to some fixed 
number, say |, is related to the saturation scale defined in the Introduction: r s — 1/Q s (y). 

We now see that there is a very close analogy between what we are describing for QCD here 
and the model that we were introducing in the previous section. So in particular, one might be 
able to formulate interaction processes in QCD with the help of a stochastic nonlinear evolution 
equation for the "gedanken" amplitude T. We already know the equation that one should get in 
the mean-field limit in which N is very large: It is the BK equation, as was rigorously proven 



above. Thus we know the equivalent of the term (u(t + At, x)) in Eq. (2.25). The noise term is 
not known, but since it is of statistical origin, it must be of the order of the square root of the 
number of dipoles normalized to N. that is to say, of order \JtJN. We may write an equation of 
the form 

9 ay T(y, k) = x(-dm k i)T(y, k) - T 2 (y, k) + a s yj2T(y, k) v(y, k), (2.29) 

where v is a noise, uncorrelated in rapidity and transverse momentum, with zero mean and unit 
variance. This equation is to be compared to the following one: 



(V, «(/..«•) = <7; </(/..»•) + ,i[t..r) - ,r(l..r) - \/ 2 u ^ x "> v(t,x), (2.30) 
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Figure 2.11: Sketch of the scattering amplitude T of a dipole of size r off a frozen partonic 
configuration. The small lines on the axis denote the dipoles ordered by their logarithmic sizes. 
Up to fluctuations, T looks like a wave front. 



which is the so-called "Reggeon field theory" equation when the noise v is exactly a normal Gaussian 
white noise, that is to say, whose non- vanishing cumulants read 



(u(t,x))=0 
(is(t, x)v(t , x )) = S(t — t)S(x — x). 



It is a stochastic extension of Eq. ( 2.23[ ). If the noise term were of the form 



2u(t,x)(l - u(t,x)) 
N 



u(t, x) 



(2.31) 



(2.32) 



instead, then this equation would be what is usually referred to as the stochastic Fisher-Kolmogorov- 
Petrovsky-Piscounov equation. The sFKPP equation and the physics that it represents is reviewed 
in Rcf. 



75 



Taking averages over events converts this equation into a hierarchy of coupled equations, which 
has a lot in common in its structure with the (modified) Balitsky hierarchy ( 2.9|2.12 ). A detailed 
study may be found in Ref. 76 . We will perform explicit calculations in this spirit within simpler 
models in Chap. [3] below. 

Based on these considerations, we may establish a dictionary between QCD and reaction- 
diffusion processes. The correspondence is summarized in Tab. |2.1| 



The mechanism for saturation of the parton densities (i.e. of the dipole number) is not known 
for sure in QCD. There is even some evidence that dipole degrees of freedom are no longer sufficient 
to describe scattering beyond some rapidity, as is understood from the appearance of sextupoles in 
the second equation of the full Balitsky hierarchy. There are also important differences between the 
reaction-diffusion model introduced above and QCD that lie in the "counting rule" of the particles 
(provided by the form of T el in the QCD case). But from the general analysis of processes described 
by equations in the universality class of the stochastic FKPP equation and the underlying evolution 
mechanisms presented in Chap. [4] we will understand that most of the observables have universal 
properties in appropriate limits, which do not depend on the details of the mechanism at work. 



We draw the reader's attention to Refs. 77 78 , where a precise stochastic equation was searched 



for in QCD. Some of the problems one may face with the use and the very interpretation of such 



equations were studied in Ref. 79 



The way in which we view high energy QCD in this review is actually not particularly original: 
It is nothing but the QCD dipole model, which was implemented numerically in the form of a Monte 
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Reaction-diffusion 
Occupation fraction u(t, x) 

Average occupation fraction (u(t,x)) 
Space variable x 



QCD 

Scattering amplitude for the probe off a frozen 
realization of the target T(y, k) 

Physical scattering amplitude A = (T) 

ln(fc 2 /A 2 ) or ln(l/a; 2 A 2 ) 



Time variable t 

Average maximum density of particles N 

Position of the front X(t) 

Branching-diffusion kernel w(— d x ) 
(oj(-d x ) = d 2 x + 1 in the FKPP case) 



ay 
1/al 

Saturation scale ln(Q 2 (y)/A 2 ) 

BFKL kernel x(-<9 lnfc2 ) 

or its equivalent in coordinate space 



Table 2.1: Dictionary between QCD and the reaction-diffusion model for the main physical quan- 
tities. A is a typical hadronic scale. 



Carlo event generator by Salam |80}|82] (see also [54] for another more recent implementation). He 
also devised and implemented a ad hoc saturation mechanism 53 that went beyond the original 
dipole model pictured in Fig. |2.6fc , but which is necessary, as we argued before. 

Before discussing more deeply the physical content of equations of the form of Eq. (2.291, we 
shall first study a model in which spatial dimensions are left out, that we will be able to formulate 
in different ways. 



2G 



Chapter 3 

The simplest saturation model 



In the previous chapter, we have understood that scattering at high energy in QCD may be viewed 
as a branching- diffusion process supplemented by a saturation mechanism. We have exhibited a 
simple toy model with these characteristics, whose dynamics is represented by an equation of the 
type fiiuftfy . 

Unfortunately, even that toy model is too difficult to solve analytically. We shall study a still 
simplified model, where there is no diffusion mechanism: Realizations are completely specified 
by the number of particles that the system contains at a given time. Of course, in this case, a 
saturation scale cannot be defined, which limits the relevance of this model for QCD. However, we 
will be able to formulate this model in many different ways, and to draw parallels with QCD. 

We start by defining precisely the model. Then, two approaches to the computation of the 
moments of the number of particles are presented. The first set of methods relies on field theory 
(Sec. 3.2). The second method relies on a statistical approach (Sec. 3.3) and will be extended in 
a phenomenological way to models with a spatial dimension in Chap. L/J We shall then draw the 
relation to a scattering-like formulation (Sec. 3.4-)- Finally (Sec. 3.5), some variants of the basic 
model are reviewed. 



Contents 

3.1 Definition [27] 

3.2 "Field theory" approach [28] 

3.2.1 Particle Fock states and their weights [25] 

3.2.2 "Pomeron" field theory [30] 
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3.4 Relation to high energy scattering and the parton model approach 1381 
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3.5.1 Allowing for multiple scatterings between pairs of particles [40] 

3.5.2 Reggeon field theory [41] 



3.1 Definition 

Let us consider a simple model in which at a given time t, the system is fully characterized by 
the number n t of particles. The evolution rules are the following. Between times t and t + dt, 
each particle has a probability dt to split in two particles. For each pair of particles, there is a 
probability dt/N that they merge into one. We may summarize these rules by the following set of 
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equations: 



n t +dt 



n t + 1 proba ntdt 

, , n t (n t -l)dt 
n t — l proba 



N 

n t proba \~n t dt 



(3.1) 



n t (n t -l)dt 
N ' 



From this, one can easily derive an equation for the time evolution of the probability P(n,t) of 
observing n particles in the system at time t: 



8P 
~dt 



(n,t) = (n- l)P(n -l,t) + 



n(n + 1) 
N 



P(n + l,t) 



n(n — 1) 
N 



P{n,t). 



(3.2) 



This is the master equation for the Markovian process under consideration. The two first terms 
with a positive sign represent the process of going from one state containing n particles to an 
adjacent one containing n + 1 or n — 1 particles respectively, while the last term simply corrects 
the probability to keep it unitary. 

By multiplying both sides of this equation by n and summing over n, we get an evolution 
equation for the average number of particles (n t ): 



d(n t ) 
dt 



= (nt) - jj(M n t - !))• 



(3.3) 



Obviously, this equation is not closed, and one would have to establish an equation for (n t (n t — 1)), 
which would involve 3-point correlators of n tl and so on, ending up with an infinite hierarchy of 
equations, exactly like in Chap. [2] for QCD (see Eqs. (2.9) and (2.12 1). 

This illustrates the difficulties one has to face before one can get an analytical expression for 
(nt), even in such a simple model. 



3.2 "Field theory" approach 

In the next subsections, we will follow different routes to get analytical results on the moments 
of the number of particles in the system at a given time t. The first one will be similar to the 
s-channel picture of QCD (see Sec. [2]), since it will consist in computing the time (equivalent to 
the rapidity in QCD) evolution of realizations of the system. The second one will be closer to the 
t-channel picture of QCD. We will see how "Pomerons" may appear in these simple systems. We 
will then examine a formulation in terms of a stochastic nonlinear partial differential equation, 
which is nothing but the sFKPP equation in which the space variable (x) has been discarded. 



3.2.1 Particle Fock states and their weights 

Statistical problems were first formulated as field theories by Doi (83] and Peliti [84j. Different 
authors have used these methods (see Ref. 185] for a review). We shall start by following the 
presentation given in Ref. (861. 



We would like to interpret the master equation (3.2 1 as a quasi-Hamiltonian evolution equation 
of the type of the ones that appear in quantum mechanics. To this aim, we need to introduce the 
basis of states |n) of fixed number n of particles. We define the ladder operators a and a) by their 
action on these states: 

a\n) = n\n - 1), a^\n) = \n + 1) (3.4) 
and which obey the commutation relation 

[a,a f ]=l. (3.5) 

The n-particle state may be constructed from the vacuum (zero-particle) state by repeated appli- 
cation of the ladder operator: 

|n> = (at)». (3.6) 
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The normalization is not standard with respect to what is usually taken in quantum mechanics. 
In particular, the orthogonal basis |n) is defined in such a way that (m\n) = n\S mn . This implies 
that the completeness relation reads 



]Tl|n)H = l. (3.7) 



We also introduce the state vector of the system at a time t as a sum over all possible Fock states 
weighted by their probabilities: 

\<t>{t))=Y,P{n,t)\n). (3.8) 



It is straightforward to see that the master equation (3.2) is then mapped to the Schrodinger-type 
equation 

^W)) = -HW)), (3.9) 

where H is the "Hamiltonian" operator 

H = {l-a})a)a-^-{l-a ] )a}a l . (3.10) 

The first term represents the splitting of particles, while the second one, proportional to 1/N, 
represents the recombination. We may rewrite H as 

H = H + H 1 , (3.11) 

where 

Ho = a^a (3.12) 

is the "free" Hamiltonian whose eigenstates are the Fock states. We now go to the interaction 
picture by introducing the time-dependent Hamiltonian 

Ui{t) = e n ° t n 1 er Uat (3.13) 

and the states \4>)i — e Uot \4>). The solution of the evolution reads 

\cf>) I =Texp(- f dt'Hi{t')] \4>o)i 

= \<t>o)i~ [ dt , H I (t , )\cf H) ) I + f dt' ( dt"H I {t')H I (t")\^ ) )i + --- 
Jo Jo Jo 

We may then compute the weights of the successive Fock states by applying this formula. Let us 
show how it works in detail by computing the state of a single particle evolved from time to time 
t, in the limit TV — oo in which there are no recombinations. We follow the usual method to deal 
with such problems in field theory. We repeatedly insert complete bases of eigenstates of Ho into 
Eq. ( |3.14| ), namely 

10)/ = |1> - f dt'Y.-^lniKniMni) + ■■■ (3.15) 



(We have kept the first two terms in Eq. (3.141 explicitely) . Using the expression for %i{t) as a 
function of Ho and Hi, together with the knowledge that the Fock states are eigenstates of Hq, 
we get 



n 1 )(n 1 \H 1 \l)+--- (3.16) 



Inserting the expression for Hi, one sees that in the infinite-TV limit, there is only one possible 
elementary transition, namely the splitting. Performing the integration over t' and computing in 
the same manner the higher orders, one finally gets the expansion 

\<t>) = e-*|l) + e-*(l - e-*)|2> + • • • + e -*(l - e"*)"- 1 ^) + • • • (3.17) 
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from which one can read the probabilities of the successive Fock states. This expansion is similar 
to the expansion in dipole Fock states introduced in Sec. [2j The n-particle states correspond to 
n-dipole states in QCD, and their weights are computed by applying successive splittings to the 



system, whose rates are given by Eq. (2.1 1. (They are just unity in the case of the zero-dimensional 
model.) 

We see that this method is well-suited to compute the probabilities of the lowest-lying Fock- 
states, and their successive corrections at finite N. But in general we are rather interested in 
averages such as (n k ), for which the weights of all Fock states are needed. We will develop a 
slightly different (but equivalent) formalism below, that will enable us to get these averages in a 
much more straightforward way. 

3.2.2 "Pomeron" field theory 

Let us introduce the generating function of the factorial moments of the distribution of the number 
of particles 

Z(z,t) =^(l + z) rl P(M)- (3.18) 



The evolution equation obeyed by Z can easily be derived from the master equation (3.2 1 



We may represent this equation in a second-quantized formalism by introducing the operators 

fe f = z, b=®-=z (3.20) 
oz 

acting on the set of states \Z) consisting in the analytic functions of z. Then we may write 

~ = -n P Z, (3.21) 

where 

n r = nl + u\, with %l = -b% u\ = -tftfb + + tf)b\ (3.22) 

A basis for the states is 

|fc) = z k , (fc| = z k (3.23) 
which is orthogonal with respect to the scalar product 

(Z 1 \Z 2 )= J d ^e-^ 2 Z 1 (z,z)Z 2 (z,z), (3.24) 

and obeys the normalization condition (A;|Z) = k\5k,i- We shall call these states "/c-Pomeron" states, 
by analogy with high-energy QCD. We may apply exactly the same formalism as before, since the 
operators b, b' have the same properties as the a, oL 

From the definition of the scalar product, it is not difficult to see that the fc-th factorial 
moment of n may be obtained by a mere contraction of the state vector \Z), computed by solving 
the Hamiltonian evolution, with a fc-Pomeron state. The following identity holds: 



(m = { ¥ ^ W J, (3-25) 

where the average in the right-hand side goes over the realizations of the system. As for the 
initial condition, starting the evolution with one particle means taking as an initial condition the 
superposition |0) + |1) of zero- and one-Pomeron states respectively. The zero-Pomeron state does 
not contribute to the evolution, hence effectively a one-Pomeron state is like a one-particle state. 
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(a) 



(b) 



(c) 



(d) 



Figure 3.1: Propagator and vertices for the Pomeron field theory. Time flows from the top to the 
bottom. 



In order to simplify the systematic computation of these moments, we may use a diagrammatic 
method and establish Feynman rules. To this aim, we write the contribution of the graphs with 



Z-vertices (corresponding to the term of order I in the expansion of Eq. ( 3. 14 1 ) , starting with a 
one-Pomeron state: 



(k\Z) 



D(-l) 1 f dt x f U dt 2 --- f'^dtt V (k\n l )~(n l \H r I \n l _ 1 )---~(n 1 \H r I \l). (3.26) 

Jo Jo Jo „~... ml ml 



Each matrix element that appears in this equation is associated to a vertex, and propagators 
connect these vertices. We read on the expression for the Hamiltonian ( |3.22[ ) that there is one 
propagator and three vertices in the theory: one splitting vertex (1 — > 2), one recombination 
(2 — >• 1) and a 2 —> 2 elastic diffusion vertex. 

The method to compute the 1 to k Pomeron transition amplitude is standard. First, one draws 
all possible diagrams for this transition that contain I vertices, including all possible permutations. 
(Note that a splitting may occur in k different ways, if k is the number of Pomerons before the 
splitting; A recombination instead may occur in k(k— 1) /2 ways). Then, the propagators (Fig. |3.1^ t) 
are replaced by 

(l\e- m °\l) = e\ (3.27) 

(where t is the time interval they span). The n-Pomeron state propagates as (n|e~ tw °|n) = e nt . 
Intermediate times are eventually integrated over. As for the vertices (Figs. [3~Tj p-d), the following 
factors have to be applied: 

(l->2): -1; (2->l): ~; (2^2): ~. (3.28) 

In addition, there is a (_i)# vertlces factor. Finally, an overall k\ factor leads to the expression for 
the factorial moment (n(n — 1) • • • (n — k + 1)). 

The lowest-order diagram for the average particle number, consisting in a simple propagator, 
reads (n) — e*. We now understand that this method leads to a more straightfoward computation 
of the moments of the number of particles than the one based on the computation of the proba- 
bilities of successive Fock states, for a single Pomeron already resums an infinity of particle Fock 
states. The Pomeron in this case is exactly like the BFKL Pomeron introduced in Sec. [2] which 



leads to an exponential increase of the scattering amplitudes with the rapidity (Eq. (3.27)). 

We now move on to the computation of higher-order diagrams in which recombinations are 
absent. First, let us recover simple results by taking the infinite- N limit. We consider the diagrams 
in Fig. |3.2| which are the only ones that survive at infinite N in the evaluation of the moment 
(n(n — 1) • • • (n — k + 1)). Using the Feynman rules, we get for each individual diagram 



t rt 

r?f„P~*2 . . . / M, O-t* = 

./,, ./,, , kl 



(-if x(-lfxe ki / dhe- 11 I dt 2 e~ t2 --- [ dt k e- tk = \-e kt (l - e"*)* \ (3.29) 
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(k! such diagrams) 

Figure 3.2: Diagrams contributing to the one Pomeron — > fe-Pomeron transition, which gives the 
moments (n(n — 1) • • • (n — k + 1)) at leading order in a 1/N expansion. 






(a) (b) (c) (d) (e) 

Figure 3.3: Diagrams up to order 1/N 2 contributing to the average of the number of particles in 
the system after an evolution over the time interval t. 

There are kl such diagrams (corresponding to all possible permutations of the Pomerons), and 
there is an extra overall k\ factor to be added in order to get the relevant factorial moment: 

(n(n - 1) • • • (n - k + 1)} = kle kt (l - e~ t ) k ~ X . (3.30) 

Next, we would like to perform the computation of the one-Pomeron — > one-Pomeron transition 
(which provides the value of (n)) within the full theory, including the recombinations. Some of 
the lowest-order diagrams are shown in Fig. |3.3| A straightforward application of the Feynman 
rules edicted above leads to the following results for the graphs that are depicted in Fig. |3.3| 



( n )ltree, Fig 


■ |33K 


= e* 
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= -2!— ( 


\-e-\l + t)) 








p 3t 
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+ 4e _t (l -t)- e 


~ 2t (2t + 5)) 






e 2t 






( n ) 2 loops, Fig 


|OJl 




-3 + e - t (4+2t- 
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( n ) I2 loops, Fig 


.[33]= 


= 4 ^<'- 


-2 + e -* {t + 2)) 





We may classify these different contributions according to their order in e t /N: We see on the 
explicit expressions that the leading terms for large t and e t /N ~ 1 are always of the form 
N(e t /N) 1+ ^ loops . It turns out that we may compute easily these dominant terms at any number 
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of loops. They stem from the graphs in which all splittings occur before all recombinations (such 
as 



3.3 d and |3.3| :). These terms build up a series that reads 

ki 



= EC" 1 )*" 1 * 1 AS=T ( 3 - 32 ) 



fc=i 

This series is factorially divergent, but is easy to resum with the help of the Borel transformation. 
Indeed, using the identity 

r+oo 

k\= dbb k e-\ (3.33) 

Jo 



replacing it in Eq. (3.321, then exchanging the integration over b and the sum over the number of 



Pomerons k, one gets 

{n)=N' 2 e- t db T e~ Ne h = N (l - Ne Ne T(0, Ne'*)) , (3.34) 

where T is the incomplete Gamma function. 

This result was obtained for the first time using a diagrammatic method in Ref. 87 . The 
authors of the latter paper also computed the next-to-leading order, that is to say, the terms of 
relative order 1/N after the resummation has been performed. The equivalent of the diffractive 
processes known in QCD were also investigated by these very authors in Ref. (88] . More results 
were obtained on that kind of models by another group in Ref. |89[|90| , using different techniques, 
which go beyond the perturbative approach. Remarkably, the latter calculations can be applied 



to some extent to QCD 91 92 . 



3.2.3 Stochastic evolution equations 

The model may also be formulated in the form of a stochastic evolution equation for the number 
of particles n t it contains at each time t. The most straightforward way of doing this would be to 



first compute the mean and variance of n t +dt given n t] with the help of the master equation (3.2 1. 
This would enable one to write the time evolution of n t in terms of a drift and of a noise of zero 
mean and normalized variance, namely: 



dn t n t {n t -l) / n t {n t -l) 

ltt =n t + ^n t + u t+dt , (3.35) 



where v is such that (vt) = and {v t v t ') — S(t — t'). This equation is similar to Eqs. (2.291 



and (2.30), except for it does not have a spatial dimension where some diffusion could take place. 
The noise term is of order y/n, as it should according to the argumentation of Chap. [2] Note that 
the distribution of v depends on n t and is not a Gaussian. This last point is easy to understand: 
The evolution of v t is intrinsically discontinuous, since it stems from a rescaling of n t , which is 
an integer at all times. A Brownian evolution (i.e. with a Gaussian noise) would necessarily be 
continuous. For completeness, we write the statistics of ft+du which is deduced from the evolution 
of n: 

[ Tdt ~ £ P roba n t dt 
vt+dt = { - £ proba 1-mdt- dt (3.36) 

l-^S-f proba ^Mdt, 

where A = n t — " t ( T ^.~ 1 ) an d a = J Uf _[_ ^("t- 1 ) ^ xhere are jumps, represented by the large 
terms proportional to l/dt. 

This formulation is not of great interest, neither for analytical calculations nor for numerical 
simulations, since it is much easier to just implement the rules that define the model in the first 



place (Eq. (3.1 1) in the form of a Monte Carlo event generator. 
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There is a better way to arrive at a stochastic evolution equation for this model, although it 
is a bit more abstract. (It is actually equivalent to the Pomeron field theory formulated before.) 
Instead of following states with a definite number of particles like above, we may introduce coherent 
states 



— z+za' 



|o>, 



(3.37) 



where z is a complex number. For real positive values of z, the state \z) is nothing but a Poisso- 
nian state, which is a superposition of |fc)-particle states, where the weight of each term follows 
the Poisson law of parameter z. For the simplicity of the argument, let us restrict ourselves to 
Poissonian states. By applying the Hamiltonian H (defined in Eq. (3.10)) to a Poissonian state 
\zt), one gets a new state \4>t+dt)'- 



dtn\z t ) 



(3.38) 



Of course, that new state is not itself a Poissonian state in general, but may be written as a 
superposition of such states. One writes 



t+dt) = I dzf(z)\z) = / dzf(z)^2e 



(3.39) 



The idea is to interpret the weight function f(z) as the probability to observe a given Poissonian 
state \z). Hence the evolution is viewed as a stochastic path 

> z t -dt z t ->• zt+dt z t+2 dt -> • • ■ (3.40) 

with well-defined transition rates from one Poissonian state to the next one. Inserting the explicit 
expression for the Hamiltonian (3.10 1 and the decomposition (3.39 1 in Eq. (3.38 1, one gets for each 
Fock state \n) 



dze~ z f(z)- 



— = e~ Zt — 

I ' TJ.I 



(n-1)! (n-2)! 



1 

N 



(n-1)! (n-2)! 



(3.41) 



Finally, this equation can be inverted for f(z) by a weighted integration over n, \ jj^z,^, , along 
an appropriate contour in the complex plane. After some straightforward algebra, we get 



f{z t +dt) = $(z t +dt - Zt) + dt ( Zt - -j ) 5'{z t+d t - z t ) + 



N 



2dt ( z t - 

1 ' N 



5"(z t+ dt - z t ) 



(3.42) 



This is a representation for the Gaussian law centered at z t + dt(z t — ^) of variance 1dt{z t — j^). 
Introducing a normal Gaussian noise v t which satisfies 

{u t ) = and (v tVt ,) = S(t - t'), (3.43) 

we may write 



dz t 
~dt 



z t 



N 



Vt+dt 



(3.44) 



where one must be careful to take the noise at time t + dt, and hence, this equation is to be 
interpreted in the Ito sense. If z t= o is a real number between and N, then the equation keeps it 
in this range. But of course the equation is valid for more general coherent states, with complex 
z t - 

This equation is suitable for numerical simulations: One may discretize the time in small steps 
At <C 1 in which case v t is distributed as 



V{vt) 



1 



V2irAt 



exp 



2At 



(3.45) 
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(In many cases, one has to use more sophisticated methods, see e.g. Ref. 193] ; A more rigorous and 
general derivation of this stochastic formulation may be obtained from a path integral formalism 
starting from the Hamiltonian (3.101, see Ref. |85|). Analytical manipulations of this equation 



using Ito's calculus are also quite easy. We are going to give an example of such a calculation 
below, avoiding unnecessary formalism. (We refer again the reader to |93 | for a textbook on a 
more mathematical handling of stochastic equations). 



We may transform the stochastic equation (3.44) to a hierarchy of equations for the factorial 



moments of the number of particles, using the relation 

(z t fe ) = (n(n-l)..-(n-fc + l)) 



= n<*>. 



(3.46) 



First, let us write Eq. (3.441 in a discretized form 



z t+ dt = zt + dt [z t 




(3.47) 



We then take the fc-th power of the left and the right-hand side, and we average the result over 
realizations. Expanding in powers of dt for small dt, we get 



(z* +dt ) = {z*)+dtk(z k t 



N 



dtk < 



.k-i 





(3.48) 



We have factorized the averages over the time intervals [t,t + dt] and [0,t], since the noise v is 
uncorrelated in time. The term proportional to dt vanishes thanks to the fact that vt+at averages 
to zero. One may think that the next term could be neglected for it is apparently proportional to 
dt 2 . Actually, it gives a contribution of order dt, be cause for discretized t, (v 2 +dt ) = l/dt. The 
dots stand for terms of order dt 2 at least. Using Eq. (3.46) to identify the factorial moments of n, 
we eventually get 



dn^ 



dt 



= k 



,0) 



N 



k(k - 1) 



,(fc-i) 



N 



(3.49) 



This equation is similar to the (modified) Balitsky hierarchy in high-energy QCD. Indeed, let us 
write explicitly the equations for the first two moments: 



d(n(n 



d(n) 
dt 
1)) 



— (n(n 



dt 



2 1 



1 

N 



1)>, 
n(n — 1)) 



N 



(n(n- l)(n - 2)) +2(n). 



(3.50) 



We note the similarity in structure with Eqs. (2.9), (2.12), except for the term 2{n) in the right- 



hand side of the second equation. This term stems precisely from the particle recombinations, and 
was absent in the B-JIMWLK/BK formalism. 



3.3 Statistical methods 

The field theory methods presented above provide a systematics to solve the evolution of the 
system to arbitrary orders in 1/N, at least theoretically. (In practice, identifying and resumming 
the relevant diagrams becomes increasingly difficult). However, it would look quite unreasonable 
to get into such an involved formalism if one were only interested in computing the lowest order in 
a large-iV expansion. Indeed, as we shall demonstrate it below, in the case of this simple model, 
an intuitive and economical calculation leads to the right answer (94] . We work it out here because 
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Figure 3.4: [From Ref. j94j] Ten different realizations of the stochastic evolution of the zero- 
dimensional model (dotted lines; N — 5 x 10 3 ). All realizations look the same, up to a shift 
in time. They are all parallel to the solution to the mean-field equation ( 3.52[ ) (dashed line). 
Note the significant difference between the latter and the average of the particle number over the 
realizations (full line). 



this line of reasoning is at the basis of the solution of more complicated models, closer to QCD, 
that we shall address in the next chapter (Chap.|4|. 

As before, we denote by n t the value of the number of particles in a given realization of the 
system. We further introduce p«(t) the distribution of the times at which the number of particles 
in the system reaches some given value n for the first time, and (n t \nt) the conditional average 
number of particles at time t given that there were ni particles in the system at time t. One may 
write the following factorization formula: 

/>oo 

(nt) = / dt P n(F)(n t \nt). (3.51) 
Jo 

This formula holds exactly for any value of h. In particular, if N is large enough, one may choose 
n such that 1 <C n <C N. 

Observing at a few realizations generated numerically (Fig. |3.4| , one sees that the curves that 
represent n t look like the solution to the mean-field equation obtained by neglecting the noise 



term in Eq. (3.35), up to a translation of the origin of times by some random to- (The curves look 
also slightly noisy around the average trend, but the noise would still be much weaker for larger 
values of N.) This suggests that once there are enough particles in the system (for nt > n 3> 1), 
the evolution becomes essentially deterministic and in that stage of the evolution, the noise can 
safely be discarded. Thus stochasticity only manifests itself in the initial stages of the evolution, 



but in a crucial way. Indeed, as one can see in Fig. 3.4 after averaging, (nt) differs significantly 
from the mean-field result, and this difference stems from rare realizations in which the particle 
number stays low for a long time. Therefore, in individual realizations, stochasticity should be 
accurately taken into account as long as n t < n. Fortunately, when the number of particles in the 
system is small compared to the parameter N that fixes the typical maximum number of particles 
in a realization, the stochastic evolution is essentially governed by a linear equation. 

Thanks to this discussion, we may assume that the evolution is linear as long as there are less 
than h particles in the system and deterministic when n t > h. It is then enough to compute Pn(t) 
for an evolution without recombinations, and (n t \n^) for an evolution without noise. The second 



quantity is most easily computed by replacing the averages of powers of nt in Eq. (3.3 1 by n 



,MF 
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and discarding the term of order 1/N. One gets a closed equation for n^ F in the form 

j MF / MF\2 

= mf _ {^t_L (3 52) 

dt N K ' 

which is solved by 

n 

where the initial condition has been chosen in such a way that n — n. 

As for the distribution Pn{t) for the waiting times t to observe n particles in the system, its 
derivation is a bit more subtle. 

Let us introduce R(n,t) the probability distribution of the first passage time at the given 
population size n, starting with n individuals at time 0. The probability Pn(t) we are looking for 
is nothing but R(l,t). 

We now establish an evolution equation for R. Recall that the evolution equation for P was 
obtained by considering the variation in the number of particles in the system between times t 
and t + dt. Here we consider the beginning of the time evolution, between times and dt. The 
probability that the system has n particles for the first time at t + dt starting with n particles at 
time 0, R(n, t + dt), is the probability ndt that the system gains a particle between times and dt 
multiplied by R(n + l,t), minus a unitarity-preserving term. In this way, after having taken the 
limit dt — >• 0, we get 

^ = n (R(n + 1, t) - R(n, t)) , with the condition R(n, t) = 8(t). (3.54) 

This equation is valid when we neglect recombination processes, which is the relevant approxi- 
mation here since we stick to the dilute regime. In order to find a solution, we introduce the 
generating function for the moments of n: 



oo 



G(«,t) = Vu n %t) (3.55) 



n=0 



and the Laplace transform 

r- + oo 



G{u,s)= dte- st G(u,t). (3.56) 
Jo 

The evolution of R implies the following equation for G: 

(1-u)^- = (s+-)g. (3.57) 



du \ M 

This equation is straightforward to integrate. Its solution reads 

G(u,s) = Cu(l - u)- 1 - 3 . (3.58) 

The constant C must be determined from the initial condition, namely from the equation R(n, t) = 
S(t), which after Laplace transform reads R(n,s) = 1. The latter means that the n-th order in 
the expansion of G in powers of u should be set to one, which writes 

% + l )r (n) = (3 ' 59) 

For large n, the Stirling formula enables one to cast the equation in the simplified form C = 
T(s + l)n~ s . We then see that R(l,s) = L(l + s)h~ s . The inverse Laplace transform of this 
function is just the Gumbel distribution: 

R(l,i) = p n (t) = ne-*- ne ~\ (3.60) 
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Plugging Eqs. (3.601 and (3.53) into Eq. (3.51), we get for the average number of particles after 

(3.61) 



t time units of evolution: 



N 



di- 



ne 



1 



N 



-(*-*) ' 



Because the Gumbel distribution is strongly damped for t < 0, the lower integration boundary 
may safely be extended to — oo. Indeed, it is easy to see that a conservative upper bound for the 
contribution of the domain ] — oo, 0] to the integral is e~ n , which is very small in the limit 



Finally, we perform the change of variable b 



ne 



/N to arrive at the form 



db- 



1 



„-ATe _ 



(3.62) 



It can be checked that it is exactly the form found through the diagrammatic approach to Pomeron 



field theory (compare Eq. (3.621 to Eq. (3.34)) 



The factorization in Eq. (3.511 and the convenient approximations that it subsequently allows 



are actually very important. Indeed, we realized that we may write the average number of particles 
at time t, whose expression would a priori be given by the solution of a nonlinear stochastic 
differential equation, by solving two much simpler problems. The key observation was the following. 
When the number of particles in the system is low compared to the maximum average number 
of particles N allowed by the reaction process, then the nonlinearity is not important, but the 
noise term is instead crucial. On the other hand, when the number of particles is large compared 
to 1, then the noise may be discarded, but the nonlinearity of the evolution equation, which 
corresponds to recombinations of particles, must be treated accurately. From this method, one 
gets an expression for (n t ) for any time up to relative corrections of order 1/N. 

When we address the problem of reaction-diffusion with one spatial dimension, we will rely on 
the very same observation. It is essentially the latter which will enable us to find analytical results 
also in that case. 



3.4 Relation to high energy scattering and the parton model 
approach 

So far, we have focussed on the factorial moments of the number n of particles in the system. 
We have seen how they may be computed from "Pomeron" diagrams, which are quite similar to 
the diagrams that appear in effective formulations of high-energy QCD. However, the relation to 
scattering amplitudes, which are the observables in QCD, may not be clear to the reader at this 
stage. In particular, we do not understand yet what would correspond to boost invariance of the 
QCD amplitudes. The aim of this section is to try and clarify these points. 

Let us consider a realization of the system of particles, evolved up to time t, that we may call 
the projectile. A convenient formalism to compute the weights of Fock states was presented in 
Sec. |3.2.1| At time t, the system of particles scatters off a target consisting of a single particle, 
and can have at most one exchange with the target, which "costs" a factor 1/N. All the particles 
in the system have an equal probability to scatter. Hence the probability that the system scatters 
reads T = n t /N. 

This way of viewing the evolution of the system makes it obviously very similar to the QCD 
dipole model introduced in Chap. [2j provided one identifies the number of particles to the number 
of dipoles and the time to the rapidity variable. The average of T over realizations is the elastic 
scattering amplitude. 

From this analogy, there is a property similar to boost invariance that should hold. Instead of 
putting all the evolution in the projectile, we may share it between the projectile and the target. 
Let us call nt> the number of particles in the projectile at the time of the interaction, and mt-t' 
the number of particles in the target. The total evolution time is the same as before. To establish 
the expression for T in this frame, it is easier to work with the probability S = 1 — T that there 
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t-t 



Figure 3.5: Representation of the scattering of two systems of particles. The systems evolve in 
time from the left to the right. The horizontal lines represent the particles, and the vertical dashed 
lines the interactions between the systems. Each of the elementary scatterings comes with a power 
of 1/N. Note the strong similarity with the QCD diagram in Fig. 2.6 1, except that in the present 
case, recombinations are included in the evolution of each of the systems. 



is no interaction. If any number of interactions were allowed between each pair of particles from 
the projectile and the target, then one would simply write S = exp(— n t rm t ^ t i /N). But since the 
number of scatterings should be limited to one per particle, one has to decrease n and m for each 
new power of 1/N, i.e. for each additional rescattering: 

S = 1 - j^nm + ^[n(n - l)][m(m - 1)] - ^^i n ( n - l)(n - 2)][m(m - l)(m - 2)] ■ • • 

(3.63) 

where the time dependences are understood, in order to help the reading. This is like a "normal 
ordering" of the expression to which we would arrive by assuming any number of exchanges. Note 
that S is not necessarily positive in a given event, and hence looses its probabilistic interpretation 
once one has performed the normal ordering. 
Taking the average over realizations, one gets 



. _v^/ n\ \ I m! \ (-l) fc 

Let us analyze this expression. 

First, if t' = t — t' ("center-of-mass frame"), the first two factors in each term of the series are of 
course identical after averaging. The sum runs over the number of actual exchanges between the 
probe and the target. A realization of the evolution, which would correspond to an event in QCD, 
is represented in Fig. |3.5[ Note that the figure is very similar to Fig. |2.6^ , except that particle 
mergings are allowed, while they have not been properly formulated in QCD yet. 

Second, this expression should be independent of t'. It is not difficult to check that this is 
indeed true by taking the derivative of (S) with respect to t' . Expressing the averages of the 
factorial moments of the number of particles with the help of the probability distributions P(n, t') 
and P(m, t — t') respectively, each term of the sum over k and m, n reads 



d(S) 



dt> 



, fixed - {PnPm PnPm) (n =1)1 (^Tfc)i ■ (3 - 65) 



The time dependence is again implicit, and we introduced the notation P n = P(n, ■) and P n = 
d t P(n, •) to get a more compact expression. The time variabe that should be used for each factor 
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is unambiguous since it is in one-to-one correspondence with the particle number index. We may 



use the master equation (3.2) to express the time derivatives: 



n „ n(n + 1) „ 
(n-l)P»_i+ v > P n+l 



n(n — 1) 
N 



Pn 



P„ 



[n-H-m]. (3.66) 



Recalling that there are sums over m, n and k which go from to oo, one may shift first the 
indices m and n in order to factorize P n P m in each term. The factors 1/N may then be absorbed 
by shifting k for the relevant terms. Then cancellations occur between the terms of both squared 
brackets in such a way that once the summations over n, m and k have been performed, the global 
result is 0. This proves the independence of (S) upon t' , that is, "boost invariance" in a quantum 
field theory language. 

We have seen that we may formulate scattering amplitudes in the zero-dimensional toy model, 
exactly in the same way as in QCD. We have seen in particular how crucial it is to include particle 
mergings consistently with the form of the interaction between the states of the projectile and of 
the target at the time of the interaction, in order to get a boost-invariant amplitude. 



3.5 Alternative models in zero dimension 

For the sake of completeness, we shall now construct some variants of the zero-dimensional model 
introduced above, since the latter were also discussed in the literature. We review two of the most 
popular models. 

3.5.1 Allowing for multiple scatterings between pairs of particles 

Instead of assuming that there is at most one single exchange between each pairs of partons, one 
may allow for any number of exchanges. Then the definition of S is modified as follows: 

(S) = (e'™) = P(nJ)P{m,t - i!)e-™ . (3.67) 



One sees immediately that if the probabilities P satisfy the master equation (3.2 1, then this 



expression cannot be boost-invariant (i.e. independent of t'). Indeed, if Eq. (3.2 1 holds, then 

P(n, t -> oo) = 5, l:N and P(n,t = 0) = <5„,i . (3.68) 
It follows that in the frame in which the projectile is at rest, 

(S}t'=o,t-K» = ( 3 - 69 ) 

while in the center-of-mass frame (if the projectile and the target share an equal fraction of the 
evolution) , 

(S) t , =ht ^ 0o = e- N (3.70) 

which is very different. Actually, in this model, the average number of particles cannot saturate 
at a fixed value N. It would not be compatible with boost invariance. 

In order to preserve boost-invariance, one has to modify the master equation. We may write 
the following general equation for the evolution of the probability: 

Pn = J2( a n-kPn-k ~ a k n P n ). (3.71) 

The coefficients a„ are the transition rates from a (n — fc)-particle state to a n-particle state. We 
determine the from the boost-invariance requirement. Actually, only the coefficient a k=1 is 
needed in the case of this model. 
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Using the same method as the one employed for checking the boost invariance in the previous 
model, we write 

d(S) 



dt' 



J2(PnPm-PnPm)(e- ! ^), (3.72) 



and express P n ,P m with the help of the master equation. Requiring that the sum over n and m 
in the right-hand side cancels leads to the rates 

\ = N (l - e-"/ N ) , (3.73) 



where the overall constant is determined from the rate in the unsaturated version of the model, 
which should hold for small values of n <C N. This model was first proposed by Mueller and 
Salam |53| . 

We see that the saturation mechanism is quite different than in the previous model. Indeed, 
the average number of particles in the system keeps growing, but at a rate that slows down and 
depends on the number of particles in the system itself. Unitarity of the scattering probability T 
is ensured first by multiple scatterings rather than by the saturation of the number of particles to 
a constant number N (up to fluctuations) . 

This model was studied in detail in Ref. (95] . The conclusions drawn in there is that the 
saturation mechanism implied by the above model is likely to be quite close to the one at work in 
QCD. We could get analytical results for this model using one of the methods presented above. In 
particular, the statistical method outlined in Sec. |3.3| would apply and lead in a straightforward 
way to the expression for (n), up to corrections of relative order 1/iV. 

3.5.2 Reggeon field theory 

Starting from the field theor y formulation in Sec. |3.2.2| we may discard the 4-Pomeron vertex 



(term (b r ) 2 b 2 /N in Eq. (3.22 1). The new Hamiltonian then reads 



H RFT = -tfb-{tf) 2 b+^b 2 . (3.74) 



The stochastic formulation reads 

dz 



d = z - — + V2z ut+dt (3-75) 



(Compare to Eq. ( 3.44 1 ) . This is the zero-dimensional version of the stochastic equation defining 
the so-called Reggeon field theory, which was intensely studied in the 70 's as a pre-QCD model 
for hadronic interactions. 

This model has peculiar properties if one insists on interpreting it as a particle model. Indeed, 



the Hamiltonian (3.741 corresponds to a generating function for the factorial moments of the 



number n of particles in the system at a given time t that satisfies 

9Z(z,t) _ dZ{z,t) z d 2 Z(z,t) 

dt - Z ^ + Z > dz N dz 2 t3.7b) 

and the corresponding master equation, obeyed by the probability P(n, t) to find n particles in 
the system at time t, writes 

dP ^ Q = -nP(n, t) + (n- l)P(n - 1, t) 
at 

+ i(n + l)(n + 2)P(n + 2,t) - — n{n + l)P(n + 1, t). (3.77) 

One can read off this equation the rates for particle creation/disappearance. One has a 1 — > 2 
splitting, with rate dt; a 2 — > annihilation with rate dt/N; and a 2 — > 1 recombination with rate 
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—dt/N. This is a negative number, and of course, it is unacceptable for a physical probability 
not to take its values between and 1. But we should not reject a priori negative probabilities 
as a formal calculation tool |96 , as long as the physical probabilities are well-defined. However, a 
Monte-Carlo code based on these negative rates turns out to be extremely unstable, and thus of 
no practical useQ 

Note that the statistical approach teaches us that in the large N ^> 1 limit, the moments of 
the number of particles in the system should not be very different than for the model with 3 and 
4-Pomeron vertices, since it is essentially the form of the fluctuations in the dilute regime which 
determines the moments at all times. 

A detailed study of the special properties of this model as well as a comparison with reaction- 
diffusion- like models may be found in Ref. [97] , 



1 We thank Al Mueller and Bo- Wen Xiao for interesting discussions on this topic, and Krzysztof Golec-Biernat 
for having brought Ref. |96| to our attention. 
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Chapter 4 

General results on stochastic 
traveling-wave equations 



In chapter^ we have shown the relevance of the stochastic FKPP equation for high-energy QCD. 
The latter represents (classical) particle models that undergo a branching- diffusion process in one 
dimension, supplemented by a saturation mechanism. Chapter^ was dedicated to a detailed study, 
from different points of view, of simplified models obtained from the former ones by switching off 
diffusion. We now go back to the study of one- dimensional models. We proceed by steps: First, we 
shall address the deterministic FKPP equation (which is equivalent to the BK equation in QCD) 
(Sec. 4.1 ). Second, we shall introduce fluctuations to get solutions for equations in the universality 
class of the sFKPP equation (Sec. ^.2 and J^.3). 



Contents 

4.1 Deterministic case: the FKPP equation 1431 

4.1.1 General analysis and wave velocity 

4.1.2 Diffusion equation with a boundary and the approach to the asymptotic 
traveling wave [46] 

4.1.3 Discrete branching diffusion [50] 

4.2 Combining saturation and discreteness 1531 

4.3 Beyond the deterministic equations: Effect of the fluctuations . . . 1571 

4.3.1 Phenomenological model and analytical results [57] 

4.3.2 Numerical simulations [60] 



4.1 Deterministic case: the FKPP equation 

We address the simplest reaction-diffusion equation, namely the FKPP equation 

d t u = d 2 x u + u-u 2 . (4.1) 

This equation was found to describe scattering in QCD under some assumptions, see Chap. [2] 

It is a mathematical theorem |98| that this equation admits traveling waves as solutions, that 
is to say, solitonic-like solutions such that 

u(t,x) =u{x-vt) (4.2) 

where v is the velocity of the wave, u is a front that smoothly connects 1 (for x — > — oo) to 
(for x — > +co). The velocities of the traveling waves and their shapes for large x are also known 
mathematically. Starting with some given initial condition which itself is not necessarily a traveling 
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wave such as Eq. ( |4.2| (but which satisfies some conditions, see below), the FKPP equation 
turns it into a stationary wave front at large times, namely a function which may be written in 
the form (4.2). The front velocity during this phase may also be predicted asymptotically. We 
informally review these results in this section. 



4.1.1 General analysis and wave velocity 



The FKPP equation (4.1l encodes a diffusion in space (through the term d x u in the right-hand 
side), a growth (term u), and a saturation of this growth (term — u 2 ). It admits two fixed-points: 
the constant functions u(t,x) = and u(t, x) = 1. A linear stability analysis shows that is 
unstable, while 1 is stable. Indeed, thanks to the growth term u in the right-hand side, a small 
perturbation u(t,x) = e <C 1 grows exponentially with time. On the other hand, a perturbation 
near 1 of the form u(t, x) = 1 — e goes back to the fixed point 1 through evolution. Hence the 
FKPP equation describes the transition from an unstable to a stable state. Therefore, we expect 
that the linear part of the equation drives the motion of the traveling wave, since the role of the 
nonlinear term is just to stabilize the fixed point u = 1. 

We shall cast the linear part of the equation into a more general form: 

d t u(t, x) = u(-d x )u(t, x), (4.3) 

where d x ) is a branching diffusion kernel. It may be an integral or differential operator. An 
appropriate kernel is, in practice, an operator such that the "phase velocity" 1^(7) = cj( 7 )/ 7 (see 
below) has a minimum in its domain of analyticity. The FKPP equation is obtained from the 
choice d x ) =9^ + 1. 

Let us follow the wave front in the vicinity of a specific value of u. To this aim, we define a 
new coordinate :ewf such that 

x = x WF + vt. (4.4) 



The solution of the linearized equation (4.3 1 writes most generally 



x ) = / 77—^0(7) exp [-7(iwf + vt) + uj(j)t] , (4.5) 
Jc 2* 71 " 

where w( 7 ) is the Mellin transform of the linear kernel uj(—d x ) (and thus 7 corresponds to —d x ), 
and defines the dispersion relation of the linearized equation. ^0(7) is the Mellin transform of 
the initial condition u(t = 0.x). Let us assume for definiteness that the initial condition is a 
function smoothly connecting 1 at x = —00 to at x — +00, with asymptotic decay of the 
form u(t = 0,x) ~ e~ loX . Then 1*0(7) nas singularities on the real negative axis, and on the 
positive axis starting from 7 = 70 and extending towards +00. Let us take a concrete example: If 
u(0,x < 0) = 1 and u(0,x > 0) = e~ loX , then 1*0(7) = 1/7 + V(7o — 7)- The integration contour 
C should go parallel to the imaginary axis in the complex 7-plane and cross the interval [0,70]. 
Each partial wave of wave number 7 has a phase velocity 

«*(7) = — , (4-6) 
7 



whose expression is found by imposing that the exponential factor in the integrand of Eq. (4.5 1 
be time-independent when the velocity v of the frame is set to v = 1*0(7). 



We are interested in the large-time behavior of u(t, x). The integrand in Eq. (4.5 1 admits a 
saddle point at a value j c of the integration variable such that 

a/( 7c ) = v, (4.7) 

that is to say, when v coincides with the group velocity of the wave packet. But the large-time 
solution is not necessarily given by the saddle point: This depends on the initial condition 1*0(7). 
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7- 7c 7* 



Figure 4.1: Front velocity as a function of its asymptotic decay rate 70 (dashed curve). It has a 
minimum at 7 = 7 C . The full line represents the actual velocity that would be selected starting 
with an initial condition decaying as e~ laX for large x. If 70 = 7- < 7c (initial condition less 
steep than j c ), then the asymptotic velocity is the phase velocity of a front which has the same 
asymptotics as the initial condition. For any 70 = 7+ > 7c- the velocity of the front is the 
minimum of the phase velocity 1^(7). 



In order to understand this point, let us work out in detail the simple example of initial condition 
quoted above. The integral has two contributions for large t: 

U (t, x) = g-7o(zwF+i>t)+"(7o)* _j_ Ke -7c(2!wF+i>t)+w(7 c )t^ {'i-S) 

up to a relative 0(1) factor k. The time invariance of u(t,x) in the frame of the wave may only 
be achieved by tuning v to one of the following two values: 

(0 v = ^, (ii) v c = ^=J{ lc ). (4.9) 

7o 7c 

In the second case, v coincides with the minimum of the phase velocity u)( , y)/ r y and in particular, 
v c < v o- The relevant value of v depends on the shape of the initial condition: 

• If 70 < 7c, i.e. the decay of the initial condition is less steep than the decay of the wave from 
the saddle-point, then one has to pick the first choice (i) for the velocity. Indeed, this is 
the only one for which the first term in Eq. (4.8 1 is time-independent, and the second term 
vanishe s at large time. Due to the fact that v c < vq. choice (ii) would make the first term 
in Eq. (4.8 1 blow up exponentially, u ~ e *to(v -v e )t _ 

• If instead 70 > 7 C , then it is the second choice (ii) that has to be made. The saddle point 
dominates, and the wave velocity at large time is independent of the initial condition. 



Figure |4T] summarizes these two cases. 

The limiting case 70 = 7 C requires a special treatment. Since it is not relevant for the physics of 
QCD traveling waves, we refer the interested reader to the review paper of Ref. [72] for a complete 
treatment also of that case. 

There exists a rigorous mathematical proof of these solutions in the case of the straight FKPP 
equation |98). These results are largely confirmed in numerical simulations for various other 
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Figure 4.2: Sketch of the shape of the front according to the large- a: behavior of the initial 
condition u(t = 0, x) ~ e~ loX . Top: 70 < j c - The asymptotic shape of the initial condition is 
conserved. The relaxation of the front is fast. Bottom: 70 > j c . The asymptotic shape of the front 
is e~ lcX , and the velocity for t = 00 is v c = ui(^y c )/^/ c . The asymptotic shape is reached over a 
distance y/i ahead of the front, and the velocity at finite time is less than the asymptotic velocity 

b y Mb- 



branching diffusion kernels, including the ones of interest for QCD (see e.g. 74 99 , and Ref. (34 



100 101 for earlier simulations of the BK equation). 



Actually, in QCD as well as in many problems in statistical physics, the initial condition is 
localized or has a finite support, and hence, its large- x decay is always very steep. Thus for the 
physical processes of interest in this review, the asymptotic front velocity, that we shall denote by 
Uno for reasons that will become clear later, reads 



(4.10) 



where the last equality defines 7 C as the value of 7 for which 1^(7) = uj (-/)/"/ is minimum. Note 
that in the context of particle physics, this result was already known from the work of Gribov, 
Levin, Ryskin [14| , and was rederived later in the framework of the BK equation (35 36 102] . 

So far, we have discussed the asymptotic velocity of the solutions to the FKPP equation as a 
function of the initial condition. When the initial condition is steep enough, then the asymptotic 
front velocity takes a fixed value which is the minimum of uj(j)/^. In the opposite case, the shape 
of the initial condition is retained (see Fig. 4.2 1. We wish to know more detailed properties of the 
wave front, such as its shape and the way its velocity approaches the asymptotic velocity. There 
are several methods to arrive at this result (which is known from rigorous mathematics, see [98]). 
At the level of principle, they all rely on a matching between a solution near the fixed point u = 1, 
and a solution of the linearized equation which holds in the tail 



4.1.2 



Diffusion equation with a boundary and the approach to the 
asymptotic traveling wave 



We now come back to the original FKPP equation (4.1 ). We have seen that the nonlinearity — u 
has the effect of taming the growth induced by the linear term u, when u gets close to 1. But 
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Figure 4.3: Shape of the solution of the branching diffusion equation (4.11 ) with a moving cutoff, 
whose position is adjusted in such a way that the maximum of u(t,x) be 1 at all times. The 
solution is represented at two different times t\ and t2, showing the soliton-like behavior of the 
solution. 



nonlinear partial differential equations are very difficult to address mathematically. It may be 
much simpler to address the linear equation 

d t u = d 2 x u + u (4.11) 

supplemented with an absorptive (moving with time) boundary condition that ensures that u(t, x) 



has a maximum value of 1 at any time. We need to work out the solution of Eq. (4.11 1 with this 
kind of boundary condition. Here, we reformulate the approach proposed in the QCD context by 
Mueller and Triantafyllopoulos [36] (see also Ref. |103| for an account of the next-to-leading order 
BFKL kernel). 



A solution to Eq. (4.11 1 with the initial condition u(t — 0, x) — S(x — xq) is given, for positive 
times, by 

«(*,*) = * exp( t -(^). (4.12) 



This solution holds if the boundary condition is at spatial infinity. 
We note that the lines x of constant u(t, x) — C are given by 

x = xct + 2t — — hit — ln(C\/47r) + terms vanishing for t — > oo. (4-13) 

(We have selected the rightmost front x > Xq)- This would be the correct expression of the position 
of the front if it were enough to solve the linearized FKPP equation, to stay around some line of 
constant amplitude closing an eye on the exponential growth behind the latter line (which would 
be tamed by the nonlinearity) . The asymptotic velocity is 2, which coincides with the critical 
velocity v c of the FKPP equation discussed above. It is corrected by a logarithmic term. We will 
see that the actual solution has the same logarithm except for the coefficient. We will be able to 
get the solution by setting an appropriate absorptive boundary which will be time dependent. We 
will proceed by steps, implementing first some fixed boundary condition in order to gain intuition 
on the form of the solution. 

So if instead of the boundary condition at infinity there is an absorptive barrier at say x = X, 
i.e. if u(t, x = X) = for any t, then a solution may be found through a linear combination of 
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the latter solution with different initial conditions, in such a way as the sum vanishes at x = X. 
This is known as the method of images. It is based on the elementary observation that any linear 
combination of Eq. (4.121 also solves Eq. (4.111. From the solution with initial condition d(x — xo), 
we subtract the solution of the same equation but with initial condition S(x — (2X — Xq)), in such 
a way that this linear combination vanishes for x = X, at any time. We get 



u x (t,x) 



(x-2X + iq)^ 



(4.14) 



At this point, let us already comment that we do not expect the solution to this problem to 
represent accurately the solution to the full FKPP equation near the boundary x ~ X since 
in that region, the details of the nonlinearity must matter. So the region of interest will be 
significantly ahead of the boundary, while the starting point xq of the evolution is at some finite 
distance of the boundary: 

x - X > 1 and x - X ~ 1. (4.15) 
One may then expand the two Gaussian terms: 



, x Q -X x- X 
u x {t,x) = -= t3/2 exp ( t 



(x - X) 2 
At 



(4.16) 



But in this equation, X does not yet depend on time. We cannot implement in a straightforward 
way a time-dependent absorptive boundary. We may get to such a solution by successive iterations: 
The main trick is to go to a frame in which the solution of the branching diffusion with a boundary 
is stationary for large times. We went through the steps of this procedure in Ref. jl]. Here we 
wish to simply argue the form of the solution from the elements we have learned so far. 

As we see in Eq. (4.16), the presence of the boundary at X requires u to vanish linearly 
at x ~ X. We expect this property to be preserved when we promote A to a function of time. 
Moreover, we know from our earlier investigations that the large- X asymptotic shape is e 



-(x-X(t)) 



From Eq. (4.16 1, we see that this shape is reached diffusively; Hence there must be a factor of the 
form e -(x-x(t)) 2 /m^ 



Putting everything together, we are lead to the ansatz 



u(t,x) = Ce~ x {x - X{t))e- (x - x{t)) 



exp 



At 



(4.17) 



where C and X are constants. 

We know that the front velocity at large time is X'(t) ~ v c — 2 and we expect a logarithmic 
correction c(t) ~ hit, so we write 

X(t)=2t-c(t). (4.18) 

Inserting Eq. (4.171 and (4.181 into Eq. ( |4.11[ ) and setting x = X(t) + a (a is a constant), we arrive 
at the equation 

c'(t) [2t{a - 1) + a 2 ] = 3a (4.19) 
which means that for large a and t, c(t) ~ | Int. Hence 



X{t) = X — X2 



X = 2t lnt- 

2 



■0(1). 



(4.20) 



The latter quantity is the position of the absorptive boundary for large times, and thus also the 
position of the front. The constant X is the position of the front in the moving frame (while X(t) 
is its position in the initial reference frame). Setting X = — 1 and C = 1, the maximum of u is 
reached at x = X(t) + 1, and is indeed equal to 1. 

For large t or in the region x — X(t) <\fi which expands with time, the Gaussian factor goes 
to 1, and we see that u(t, x) only depends on one single variable x — X(t). This was expected: It is 
precisely the defining property of traveling waves. But in addition to these asymptotic solutions, 
we get from this calculation the first finite-t correction to the front shape and front velocity. 
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Actually, the speed of the front is intimately related to its shape. At time t, the front has 
reached its asymptotic shape over the distance \fi from the saturation point. This remark will be 
important in the following. 

We have derived the solution of a problem that was not exactly the initial one, however, we 
believe that the shape of the front in its forward part (u <C 1, namely for x — X(t) ^> 1) as 
well as its velocity are quite universal. Heuristically, these properties arc completely derived from 
the linear part of the equation. For this reason, the front is said to be "pulled" by its tail. The 
nonlincarity only tames the growth of u near u ~ 1, and so its precise form should not influence the 
front position itself, at least at large enough times. Thus we expect these solutions to have a broad 
validity, only depending on the diffusion kernel, and so, may be obtainable from our calculation 
up to the replacement of the relevant parameters. 



For the more general branching diffusion kernel in Eq. (4.3), the velocity of the front would 
read 

dX{t) 
dt 



7c 



3 

27ct 



(4.21) 



where j c solves u>(j c ) = 7 c w'(7 c ), as was explained in Sec. 4.1.1 The front shape in its forward 
part x — X{t) 3> 1 is represented by the equation 



u(t,x) cx (x - X{t))e-^ {x - x(t)) 



exp 



(4.22) 



up to an overall constant. Fig. |4. 3 [ represents a sketch of the solution at two different times. The 
large-time shape is an exponential decay, 



u(t, x) 



- 7c (x-X(t)) 



(4.23) 



up to a linear growth, and from Eq. (4.221, this shape extends over a range 



L = x — X(t) ~ v /2w"(7 c )t. (4.24) 
In other words, the time needed for the front to reach its asymptotic shape over a range L reads 



L 2 
2^" (7c)' 



(4.25) 



Through our simple arguments and calculation, we got the lowest order in an expansion of the 
front shape and position at large times. The next corrections to X(t) would be of order 1 (this 
constant depends on the way we define the position of the front), followed by an algebraic series in 
t whose te rms all vanish at large t. The first next-to-leading term in the series has been computed 
It turns out to be of order 1 / \ft. We will not reproduce the calculations that 



(see Ref. 



104 



lead to it because they are rather technical and there is already a comprehensive review paper 
available on the topic 72 . But let us write the result for the position and the shape of the front 
at that level of accuracy, for the more general branching diffusion kernel given by Eq. (4.3 1. To 



that accuracy, the front position reads 104 105 



X(t) 



w(7c). 



7c 




0(l/t). 



(4.26) 



For the simple FKPP case, we recall that oj(—d x ) = d 2 + 1, then j c = 1 and w(7 c ) = 2. The first 
two terms in the last equations match the ones found in Eq. (4.21 ). The shape of the front in its 
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forward part has the following form 1 104 105 : 



i(t, x) = Cie-^- x(t)) exp (-z 2 ) x 



-A-(/)) - -- - ( 3 - 2C 2 + 7 ^, (7 ( ] c) ) 2 



2 7c ^ 3 )( 7c ) ( 1 



3 w"(7c) 



+ -2F 2 [l,l;f,3;z 2 ] j : 



+ 6v^F (1 - iFi [-1, § ; z j) « + 0(1/Vt) , (4.27) 



where 

x-X(t) 

C\ and C2 are constants, and 2F2, \F\ are generalized hypergeometric functions 
the first and second lines match with the result of our calculation (Eq 



(4.28) 



The terms in 
These expressions 



g-22)) 

should apply also to QCD, up to the relevant replacements given in Tab. |2.l| 

So far, we have considered equations of the type of Eq. (4.1l as saturation equations, in the 
sense that they describe the diffusive growth of a continuous function u until it is tamed for u ~ 1. 
We will see below that these equations may actually be given a different physical interpretation. 



Relevance of this formalism to the BK equation 

In order to check that the formalism used to arrive at a solution to FKPP-like equation applies to 



the BK equation in QCD, we performed in Ref. 74 a numerical simulation of the BK equation. 
We compared the velocity of the obtained traveling wave with Eq. (4.261, using either the full 
expression with three terms or a truncation of it keeping only the two dominant terms. We 
defined the saturation scale by the equation A(y, Q s (y)) — K, and we chose different values of k. 
We see in the plot of Fig. |4.4| that the numerical result is consistent with the analytical expectations 
(using the dictionary in Tab. 2.1|, although the complete discussion is quite subtle. All details of 



the numerics and the discussion of the results were published in Ref. 1 74 . 



4.1.3 Discrete branching diffusion 

We have investigated the solutions of the FKPP equation in a mathematical way, without dis- 
cussing the physics that may lead to such an equation. The absorptive boundary that we have put 
replaces the nonlinear term in the FKPP equation, whose role is to make sure that u never exceeds 
the limit u — 1. Hence we have thought of this boundary as a way to enforce the saturation of 



some density of particles. Actually, the FKPP equation (4.1 1 may stem from a branching diffusion 



process in which the number of particles is unlimited, and thus, for which there is no saturation at 
all. As a matter of fact, this is what the BK equation describes in QCD: An exponentially growing 
number of dipolcs, stemming from the rapidity evolution of a hadronic probe, scatters off some 
target. The overall interaction probability is unitary because multiple scatterings are allowed (the 
interaction probability of n dipoles is actually of the form 1 — e _Q «"), but not because there is a 
saturation of the number of dipoles in the wavefunction of the probe. We refer the reader back to 
Fig. |2.4| for a picture of the process. 

To illustrate how the FKPP equation arises in such a simple model of branching diffusion, let 
us consider a set of particles on a line, each of them being indexed by a continuous variable x. 
(Such a model was considered for instance in Ref. [71] )• We let the system evolve according to 
the following rules. During the time interval dt, each particle has a probability dt to split in 2 
particles. Unless it splits, it moves of the small random amount 8x, which is a Gaussian variable 
distributed like 

1 ( {bxf 



p(Sx) 



V^ndt 



exp 



Adt 



(4.29) 
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Figure 4.4: Velocity of the QCD traveling wave as a function of the rapidity. The different curves 
correspond to the numerical simulation, the saturation scale being defined in various ways with 
the help of the parameter k (see the text for the definition of the latter), and to the analytical 



formula of Eq. (4.26), either truncated after the second term ("Analytic 2 terms"), or complete 



('Analytic 3 terms") (see the dictonary in Tab. 2.1 for the notations) 




Figure 4.5: Example of branching diffusion process on a line (see the text for a mathematical 
description of the evolution rules) . If the number of individuals is limited by a selection process 
which, at each new branching, eliminates the individual sitting at the smallest x as soon as the 
total number of individuals reaches say N (N = 10 in this figure), then only the branches drawn 
in thick line survive. 
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Let us consider the number of particles n(t, x) contained in an interval of given size Ax centered 
around the coordinate x. At time t = 0, the system is supposed to consist in a single particle 
sitting at the origin x = 0. A sketch of a realization of this model is shown in Fig. |4.5| From the 
evolution rules, we easily get an equation for the average number of particles (n): 



(n(t + dt,x)) =dt2(ri) + (l-dt) / d(Sx)p(Sx)(n(t,x - Sx)) 



which reads, after replacing p by Eq. (4.291 and after the limit dt — >• has been taken, 



d(n) 



d 2 (n) 
dx 2 



(4.30) 



(4.31) 



All the dependence on the size Ax of the "bin" is contained in the initial condition. It is clear that 



for large enough times, the solution to this equation is given by Eq. (4.12) 
Let us now define 

S(t, x 



_ e -n(t,x)/N 



(4.32) 



where N is some (large) constant. This definition is reminiscent of the S'-function, related to the 
scattering amplitude, introduced in the discussion of the BK equation in Chap. [2] At a fixed 
time and for large enough x, n(t,x) <C N and thus 1 — S(t,x) ~ n(t,x)/N — > 0. For any x, the 
exponential makes sure that S ranges between and 1. Thus S (or 1 — S) has the shape of a 
traveling wave. Its position X(t) is the value of x for which n(t,x) is some given constant say of 
the order of N. In the mean-field limit in which n is replaced by its average (n), it is very easy to 



compute X(t) from the form of the solution (4.121. We get (see Eq. (4.13 1) 



X(t) =2t- ~lnt 



(4.33) 



up to a constant. 

On the other hand however, the average of S over events, namely A = 1 — (S) obeys the FKPP 
equation. Indeed 



(4.34) 



(S(t + dt,x)) = dt(S(t,x)) 2 + {l-dt) J d(6x)p(Sx)(S(t,x-6x)) 
In the limit dt — > and rewriting the equation with the help of A, we get 



dA 
~dt 



d 2 A 

dx 2 



A -A 2 



(4.35) 



Hence A is a traveling wave at large times, and its position X(t) is given by Eq. (4.20). It is 
obviously behind by a term Int with respect to the value of x for which the average number 
of particles has a given constant value (compare Eq. (4.201 and Eq. (4.13 1). Furthermore, the 



probability distribution of the position of the rightmost particle (or of the fc-th rightmost particle 
for any given k) may also be derived from the FKPP equation. (Brunet and Derrida have recently 
given an advanced discussion of the statistics of the position of these particles, see Ref. 106[|107| ). 
It turns out that in any event, the average x for which n(t, x) has a given value, say no, moves with 
the FKPP velocity which can be read off from Eq. (4.201. This is much slower than the rate of 



change of X{t) when the latter is defined as the implicit solution of the equation (n(t, X(t))) = uq. 

All this may seem a bit paradoxical. But actually, it is just related to the fact that (e~ n / N ) 
cannot be approximated by e~^/ N . We may understand it in the following way. By taking the 
average of n, we have somewhat forgotten a fundamental property of n: its discreteness. Indeed, 
it only takes integer values, and in particular, the distribution of n in a realization has a finite 
support: At any time, there is a value of x to the right of which there are no particles at all. n 
obeys a stochastic equation. This is not the case for (n), which just obeys an ordinary branching 
diffusion equation. 
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<n> 




Figure 4.6: Solution of the branching diffusion equation (4.361 with a moving absorptive boundary 
that forces (n) to vanish at the point X(t) such that (n(t, X(t) — 1)) = 1. Two different times are 
represented. 



In order to recover the effect of the discreteness of n and compute the velocity, we may again 
use the absorptive boundary trick. Let us solve the linear equation 

d t {n) = d 2 x (n) + (n) (4.36) 

with an absorptive boundary. The latter will be placed in such a way that at a distance of order 



one to its left (we will focus on the right-moving wave), (n) = 1 (see Fig. 4.6 1. There is no 
difference in principle with the boundary calculation that we have performed before, except that 
the absorptive boundary is now placed to the right of the front (i.e. xq < X in the notations 
used above). Thus we find without any further calculation that the realizations of n move, on the 



average, with the FKPP velocity (4.21 1 



4.2 Combining saturation and discreteness 

We have seen that physically, the FKPP equation (or the BK equation in QCD) may be interpreted 
either as an equation for the growth, diffusion and saturation of a continuous function, or as the 
evolution equation for the average of a bounded function of a discrete (thus stochastic) branching 
diffusion process. For each of these interpretations, we may find the main features of the solutions 
by imposing one absorptive boundary on the linear partial differential equation encoding branching 
diffusion. In one case, the boundary is a cutoff that prevents u to be larger than 1: It represents 
saturation, i.e. the explicit nonlinearity present in the FKPP equation. In the other case, the 
boundary forces the function n that represents the number of particles to vanish quickly when n 
becomes less than 1 . Formally, it actually models the stochasticity due to the intrinsic discreteness 
of the number n of particles, and avoids to address a stochastic equation directly. 

In physical cases such as reaction-diffusion processes for finite N, we define u(t, x) as the 
number of particles per site (or per bin) in x normalized to N. Hence it takes discrete values: 
1/N, 2/N etc... While for large N discreteness is unlikely to play a role in the region u ~ 1, it 
is expected to be crucial when u ~ 1/N. It is thus natural to impose the two boundaries: one 
representing saturation of the particle number, the other one taking care of the discreteness of 
the same quantity. A model that these two cutoffs may represent is for example, the branching 
diffusion model in Sec. |4.1.3[ but in which the total number of particles is limited to N by keeping 
only the N rightmost ones at each new branching. It is clear that the function U(t,x) defined to 
be the number of particles to the right of some position x normalized to the maximum number N 



is, for large enough times, a front connecting 1 (for x — > — oo) to (for x — > +oo) (see Fig. 4.7 1 
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fraction of particles to the right of x 




Figure 4.7: Branching diffusion model of Sec. |4.I.3| with selection that limits the total number of 
particles to N. The function U(t, x), which is the number of particles to the right of x normalized 
to the maximum number of particles N, is represented. One sees that the fraction of particles to 
the right of x looks like a traveling wave front. 



Reaction-diffusion problems (described by nonlinear stochastic partial differential equations) 
were interpreted as branching diffusion problems taking place between two absorptive boundaries 
for the first time by Brunet and Derrida in Ref. |I08 and later, independently, by Mueller and 



Shoshi in the case of QCD in Ref. 40 



Note however that in the context of the QCD parton 

Mueller and Shoshi 



41 



model, the present interpretation of the cutoffs was only found in Ref. 
introduced the both cutoffs for reasons related to the boost-invariance of the QCD amplitude 
(The discreteness cutoff was thought as the symmetric of the saturation cutoff under some boost). 
The duality of the two boundaries, that is to say of the dense and dilute regimes of the traveling 
wave, was studied more deeply in Refs. [i09HIf3 



Before moving on to the technical derivation of the shape and position of the front in this case, 
let us figure out what we expect to find. 

Starting from the initial condition which we assume to be one or a few particles, the front 



builds up and its velocity increases with t (see Eq. (4.21 1) until it reaches its asymptotic shape, 
which is a decreasing exponential e~ lc ^ x ~ x ^ that holds for all x — X(t) 1. (X(t) is here the 
position of the bulk of the front, say for example of the leftmost surviving particle). But if the 
front is made of discrete particles, then it has a finite support, and the exponential shape may not 
extend to infinity to the right, since u(t,x) has to be either larger than 1/A, or zero. It cannot 
take values that would be a fraction of l/N in realizations, and thus, we cannot accommodate 
the shape e~ lclyX ~ x ^ for arbitrarily large values of x, since it would mean authorizing arbitrarily 
small positive values of u(t, x). From Eq. ( 4.25[ ) and from the shape of the asymptotic front (4.23 1, 
the exponential shape sets down to u = l/N at time 



^rclax 



2w"( 7c ) V 7, 



In A 



(4.37) 



Beyond, the front cannot develop any longer, and thus, its shape and velocity remain fixed. t Te \ ax 
is the time that is needed for the front to relax from any perturbation, which is why we have put 
the subscript "relax". 

From Eq. (4.211 evaluated at t = £ ro iaxj we get the new asymptotic velocity, which takes into 



account the effects of discreteness, in the form 

dX(t) _uj{ 1c ) 3 7c w"( 7c ) 



dt 



7c 



c hi A 



(4.38) 
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Figure 4.8: Sketch of the solution to the branching diffusion equation with two boundaries. 



The calculation of the constant c requires a proper account of the exact shape of the front. We 
shall turn to this calculation now. 

As announced, we are now going to solve the linear branching diffusion equation with two 
absorptive boundaries: one representing saturation, the other one discreteness. Using the intuition 
gained from the study of the deterministic FKPP equation, we write the ansatz 



(4.39) 



L is a constant which will represent the size of the front, which is essentially equal to Lq = lniV/7 c 
for large N. When u> is expanded to second order around the eigenvalue j c , then ip obeys the 
partial differential equation 

W = \®> + ^L( w '( 7c ) _ X '(t))^ (4.40) 

where we have defined 

^"{lc)t , x-X(t) . A 

V = L a and p= j±L. (4.41) 

We have only kept the dominant terms for large L. We see that ui' (j c ) — X' (t) has to scale like 1/L 2 
for all terms of this equation to be relevant, as was already guessed heuristically. The coefficient 
of 1/L 2 must be chosen in such a way that in the large-y limit, there is a nontrivial stationary 
solution. We will check that the correct ansatz is 

*'(*) = ~ + o(l/L 2 ) (4.42) 



Equation (4.40) then becomes 

d v i, = \d 2 p ^ + (4-43) 

up to higher-order terms when L is large. 

We now implement the absorptive boundaries at p = and another one at p — 1 (which 

corresponds to a distance L between the boundaries in x-coordinates, i.e. to the natural size of 
the stationary front). The boundary conditions formally read 

%jj(y, p = 0)=0 and tp(y, p = 1) = 0. (4.44) 
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As for the initial condition, for reasons that will become clear later, put a localized mass close 
to the rightmost boundary, namely we write 



^(y = 0,p) = 8(p- 1 + d) 



2 7c<5 

IX' 



(4.45) 



where a is a constant of order 1/L, and therefore a <C 1. The value of the weight e 1 " 6 actually 
corresponds to putting one single particle at a distance 5 to the right of the rightmost boundary. 
As we will see, such a weight corresponds to a fluctuation added to a stationary solution. But in 
this section, we shall only focus on the large-time behavior: The initial condition will be forgotten 
through the time evolution. 



The solution of Eq. (4.431 with the conditions (4.441 and (4.451 reads 

2e* 4 ~ 



ips(y,p) 



L 2 



y^(-l) n+1 sinTma sin 



irnpe 



(4.46) 



While the full solution with all harmonics will be of interest later, we shall discuss here only the 
stationary solution. We see that for large y, the higher harmonics are suppressed exponentially 
with respect to the fundamental mode n — 1, which gives the following contribution: 



2e 7c«5o 



— =7; — sm 7ra sm np ~ 
L 2 a«i 



27rae 7 <= <5 ° 
IX 



sm ir p. 



(4.47) 



Thanks to the choice (4.42) for X'(t), this solution has no y dependence, and leads to a stationary 
u in the frame of the front. The expression (4.471 is independent of the initial condition except 
for the overall normalization. The value of 8q, which characterizes the initial condition, will be 
adjusted later. Undoing the changes of variables which trade u for ip, x for p and t for y (Eq. (4.39 1), 
the stationary solution us reads 



ug (t, x) = e 



_ - 7c ( x -X(i))- 



2nae< Jo [ . n(x-X(t)) 



L 2 



L sin 



L 



(4.48) 



We further require that uo(t, x) ~ 1 for x = X(t) + aL, where aL is a constant of order 1. This 
condition is satisfied if we set 5q ~ 31nL/7 c . Indeed, with this choice, 



u So (t,X(t) + aL) ~ 2ir 2 aaL 



2- lc aL 



(4.49) 



Since aL and aL are constants, the right-hand side is just a number of order 1. 
All in all, the final solution reads 

u(t,x) oc K e-^ x - x ^L S m ^ X ~ X{t)) 



(see Fig. 4.8) where the size of the front is 



L 



In A 

7c 



and its velocity reads, from Eq. (4.421, 



dX(t) 
dt 



7T 2 W"(7 C ) (jj(j c ) 7T 2 7 c w"(7 c ) 



2 7c L 2 



7c 



2 In N 



(4.50) 



(4.51) 



(4.52) 



The subscript BD stands for "Brunet-Derrida" after the first authors who wrote down such an 
expression. In the FKPP case, namely for 1^(7) = 7 2 + 1, 7 C = 1 and uj("{ c ) — oj"(j c ) = 2. 
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Figure 4.9: Evolution of the front with a forward fluctuation. At time to, the primary front 
extends over a size L and is a solution of the branching diffusion equation with two appropriate 
boundaries. An extra particle is stochastically generated at a distance 5 with respect to the tip of 
the primary front. At a later time, the latter grows detcrministically into a secondary front that 
is a bit slower, and that will add up to the primary one. The overall effect, after relaxation, is a 
shift to the right of the distance R(S) with respect to the position of the front if a fluctuation had 
not occured. 

4.3 Beyond the deterministic equations: Effect of the fluc- 
tuations 

So far, we have actually solved deterministic equations although we were addressing a model with 
a discrete number of particles, that therefore had necessarily fluctuations. Our procedure gave 
the leading effects. We shall now incorporate more fluctuation effects, in a phenomenological way. 



(We shall essentially review Ref. 114 ). 



4.3.1 Phenomenological model and analytical results 

The two-boundary procedure has led to the following result: The front propagates at a velocity 



vbd in Eq. (4.52) lower than the velocity predicted by the mean-field equation (4.211, and its 
shape is the decreasing exponential e~ 7c ^ x ~ x ^ down to the position 

, . In A . . 

xtm(t) = VBvt -I , (4.53) 

at which it is sharply cut off by an absorptive boundary. This boundary was meant to make 
the front vanish typically over one unit in x, hence to implement discreteness on a deterministic 
equation. 

But since the evolution is not deterministic, it may happen that a few extra particles are sent 



stochastically ahead of the tip of the front (See Fig. 4.9 1. Their evolution would pull the front 
forward. To model this effect, we assume that the probability per unit time that there be a particle 
sent at a distance 5 ahead of the tip simply continues the asymptotic shape of the front, that is 
to say, the distribution of 6 is 

p{6) = Cxe-**, (4.54) 

where C\ is a constant. Heuristic arguments to support this assumption were presented in 
Ref. |114 . Note that while the exponential shape is quite natural since it is the continuation 
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Figure 4.10: ipSoiViP) + ^siUiP) [ see Eqs. (4.47), (4.46)] for different values of the reduced time 
variable y after a fluctuation of size 6 = 5 has occurred at y = 0. In this plot, the size of the front 
is L = 10, and a — 0.1. We see how the fluctuation, initially localized at the tip of the front, gets 
smeared uniformly over the width of the front as y gets large. Eventually, a small forward shift 
X — > X + B would be needed in order to absorb it and recover the stationary front. 



of the deterministic solution (4.221 in the linear regime, the fact that C\ need to be strictly 



constant (and cannot be a slowly varying function of ft) is a priori more difficult to argue. 

Once a particle has been produced at position x t i P + S, say at time t , it starts to multiply (see 



Fig. 4.9 1 and it eventually develops its own front (after a time i re iax of the order of L ), that will 



add up to the deterministic primary front made of the evolution of the bulk of the particles. 
Note that the philosophy of our phenomenological approach to the treatment of the fluctuations 



is identical to the spirit of the statistical approach in Sec. 3.3 developped for the zero-dimensional 
model. Whenever the number of particles is larger than n (n = 1 here), we apply a deterministic 
nonlinear evolution. Fluctuations instead are produced with a probability which stems from a 
linear equation. The difficulty here is that we are unable to solve either of these equations for 
arbitrary initial conditions and thus we have to make conjectures on the form of their solutions. 

Let us estimate the shift in the position of the front induced by these extra forward particles. 
The solution of the diffusion equation is the superposition of the large-time stationary solution 
us given by Eq. (4.481 with <5 — 31 n ^/7c- an d of the solution ug of the diffusion equation with 



the generic initial condition characterized by 8 (see Eq. (4.46 1), up to a multiplicative constant 



C-2, of order 1 that we do not control in this calculation, since it certainly depends on the detailed 
shape of the fluctuations. We write 



t (t, x) = us.it, x) + C 2 u 5 (t, x) = e -7.(*-*(*)>L [^ (y, p) + C 2 Mv, P)\ 



(4.55) 



up to the replacement of the variables by their expressions ( 4.41| . The presence of the second term 
alters the shape of the front (the front eventually relaxes back to the sine shape in Eq. (4.48 1), 
see Fig. |4.10| But of course, we want to keep the normalization condition for u, namely for 
some appropriate value of x, u is required to equate Eq. (4.49) at all t. This is possible by 



shifting the value of x at which we enforce the normalization condition from x = X{t) + aL to say 
x = X(t)+aL+R(t, S). This is equivalent to shifting the position of the front X(t) — > X(t)+R(t, 5). 
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Equation (4.55) then leads to 

u(t, X(t) + aL + R{t, 5)) = 2TT 2 aaL 2 e^ aL -^ R ^ s) 



l + C 2 



2n 2 aa L 



Equating the right-hand sides of Eq. (4.56) and Eq. (4.491, we get 



R(t, 5) = — In 

7c 



i + c 2 \ n _ - ' 



2n 2 aa L 



(4.56) 



(4.57) 



where only the lowest orders in a, a in the expansion of ifis must be kept. With the help of 



Eq. (4.461, it is then straightforward to arrive at an explicit expression of R. 



In the large time limit in which only the fundamental mode survives in the expression of ip, 
we get the shift 

R(8) = — In [l + Ci—r- . (4.58) 



7c 



The probability distribution (4.541 and the front shift (4.58) due to a forward fluctuation define 



an effective theory for the evolution of the position of the front X(t): 



X(t + dt) = 




v BD dt 
v B i)dt 



R(S) 



proba. 1 — dt J °° dSp(5) 
proba. p{5)d8dt. 



(4.59) 



From these rules, we may compute all cumulants of X(t), by writing the evolution of their gener- 
ating function, deduced from the effective theory (4.591: 



9 1 



,\X(t) 



Aw B D 



dSp(S) (t 



\R(8) 



1 



(4.60) 



The left hand-side is a power series in A whose coefficients are the time derivatives of the cumulants 
of X(t). Identifying the powers of A in the left and right-hand sides, we get 



V - v BB = 

[n-th cumulant] 

t ~ 



d6 P (5)R(S) = 
dSp(S)[R(S)] n 



CiC 2 31nL 

7c 7c£ 3 
= dC 2 n\C,{n) 

7c lc L3 



(4.61) 



We see that the statistics of the position of the front still depend on the product C\C 2 of the 
undetermined constants C\ and C 2 - We need a further assumption to fix its value. 

We go back to the expression for the correction to the mean-field front velocity, given in 



Eq. (4.52). From the expressions of R(S) (Eq. ( 4.58 1) and of p(S) (Eq. (4.54 1) , we see that the 
integrand defining V — «bd in Eq. (4.61 ) is almost a constant function of S for S < Sq = 3lnL/"f c , 
and is decaying exponentially for 5 > 5q. Furthermore, R(Sq) is of order 1, which means that 
when a fluctuation is sent out at a distance S ~ Sq ahead of the tip of the front, it evolves into 
a front that matches in position the deterministic primary front. We also notice that when a 
fluctuation has 5 < S Q , its evolution is completely linear until it is incorporated to the primary 
front, whereas fluctuations with 6 > Sq evolve nonlinear ly but at the same time have a very 
suppressed probability. We are thus led to the natural conjecture that the average front velocity 



is given by vbd in Eq. (4.521, with the replacement 



InTV . 
L -> L e g = h d 



InA lnlnTV 
+3 , 



(4.62) 



namely 



V 



7rV'( 7c ) 



27c (fe 



N i 3 In In TV 



(4.63) 
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The large-TV expansion of the new expression of the velocity yields a correction of the order of 
In In TV/ In 3 TV to the Brunet-Derrida result, more precisely 



w( 7c ) 7t 2 7 c cj''(7c) . 2 „, s31nlniV 

V = — h 7T 7 C W (7 C j - 



7c 



2 In 2 TV 



In 3 TV 



Eqs. (4.611 and (4.641 match for the choice 

dCi = 7rV'( 7 c). 



(4.64) 



(4.65) 



From this determination of C1C2, we also get the full expression of the cumulants of the position 
of the front: 

[n-th cumulant] 22/// ra!C(n) 



(4.66) 



' 7™ In 3 TV ' 

We note that all cumulants are of order unity for t ~ In 3 TV, which is the sign that the distribution 
of the front position is far from being a trivial Gaussian. This makes it particularly interesting. 
On the other hand, the cumulants are proportional to k = t/\n 3 N, which is the sign that the 
position of the front is the result of the sum of n independent random variables, and as such, 
becomes Gaussian when k is very large. The properties of the statistics of the front position were 
investigated in some more details in Ref. 1 1 1 5 . 



Thanks to our discussion in Chap. [2] we see that these results should apply to QCD with the 
relevant substitution of the kernel u; and of the parameter TV according to Tab. |2.1| 



4.3.2 Numerical simulations 

The results obtained so far rely on a number of conjectures that no-one has been able to prove so 
far. In order to check our results, let us consider again the model introduced in Sec. |2.2.2| The 
first step to take before being able to apply our results to this particular model is to extract from 
the linear part of Eq. (2.261 the corresponding function a; (7), and then to compute j c . Setting 
Ax = At = 1, we get 

w(7) = In [1 + A + Pl (e-T - 1) + Pr {e< - 1)] , (4.67) 

and 7 C is defined by uj(~f c ) = 7 c u/(7 c ). 

For the purpose of our numerical study, we set 



Pi = Pr = 0.1 and A = 0.2 . 



(4.68) 



Simulated realizations for this set of parameters arc shown in Fig. |4.11| 



From (4.67), this choice leads to 



7c = 1.352 - •• , u/( 7c ) = 0.2553- 
w "( 7c ) = 0.2267- •• . 



(4.69) 



Predictions for all cumulants of the position of the front are obtained by replacing the values of 



these parameters in Eqs. (4.64 1,(4.66 1 



Technically, in order to be able to go to very large values of TV, we replace the full stochastic 
model by its deterministic mean field approximation u — > (u) , where the evolution of (it) is given 
by Eq. (2.26 I, in all bins in which the number of particles is larger than 10 3 (that is, in the bulk of 



the front). Whenever the number of particles is smaller, we use the full stochastic evolution (2.241 



We add an appropriate boundary condition on the interface between the bins described by the 
deterministic equation and the bins described by the stochastic equation so that the flux of particles 
is conserved [116 . This version of the model will be called "model I". Eventually, we shall use 
the mean field approximation everywhere except in the rightmost bin (model II): at each time 
step, a new bin is filled immediately on the right of the rightmost nonempty site with a number 
of particles given by a Poisson law of average 6 = N(u(x, t+l)\{u(x, £)}). We checked numerically 
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Figure 4.11: 1000 realizations of the model introduced in Sec. 2.2.2 at two different times (dotted 
lines), and the average of u over the realizations (full line). One clearly sees that (u) does not 
keep its shape upon time evolution, which shows that the traveling wave property of the FKPP 
equation is lost due to the stochasticity. This point is addressed in some detail in Chap. [6] 



that this last approximation gives indistinguishable results from those obtained within model I as 
far as the statistics of the position of the front is concerned. 
We define the position of the front at time t by 



X, 



oo 



it(x, t). 



(4.70) 



x=t) 



We start at time t = from the initial condition u(x, 0) = 1 for x < and u(x, 0) = for x > 0. 
We evolve it up to time t = In 2 N to get rid of subasymptotic effects related to the building of 
the asymptotic shape of the front, and we measure the mean velocity between times In 2 N and 
16 x In 2 N. For model I (many stochastic bins), we average the results over 10 4 such realizations. 
For model II (only one stochastic bin), we generate 10 5 such realizations for N < 10 50 and 10 4 
realizations for N > 10 50 . In all our simulations, models I and II give numerically indistinguishable 
results for the values of N where both models were simulated, as can be seen on the figures (results 
for model I are represented by a circle and for model II by a cross) . 



Our numerical data for the cumulants is shown in Fig. 4.12 together with the analytical predic- 



tions obtained from (4. 64), (4. 66 I (dotted lines in the figure). We see that the numerical simulations 
get very close to the analytical predictions at large TV. However, higher-order corrections are pre- 
sumably still important for the lowest values of N displayed in the figure. 

We try to account for these corrections by replacing the factor (In N) /-f c — Lq in the denomi- 



nator of the expression for the cumulants in Eqs. (4. 64), (4.661 by the ansatz 



Jeff 



La 



31n(lniV) 

7c 



ln(lniV) 
lnA^ 



(4.71) 



The two first terms in the r.h.s. are suggested by our model. We have added two subleading terms 
which go beyond our theory: a constant term, and a term that vanishes at large N. The latter 
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Figure 4.12: [From Ref. |114| ] From top to bottom, the correction to the velocity given by 
the cutoff theory and the cumulants of orders 2 to 5 of the position of the front in the stochastic 



model. The numerical data are compared to our parameter- free analytical predictions (4. 64), (4. 66 1, 
represented by the dashed line. 



are naturally expected to be among the next terms in the asymptotic expansion for large N. We 
include them in this numerical analysis because in the range of N in which we are able to perform 
our numerical simulations, they may still bring a significant contribution. 

We fit (4.71 ) to the numerical data obtained in the framework of model II, restricting ourselves 
to values of N larger than 10 30 . In the fit, each data point is weighted by the statistical dispersion 
of its value in our sample of data. We obtain a determination of the values of the free parameters 
c = -4.26 ± 0.01 and d = 5.12 ± 0.27, with a good quality of the fit ( X 2 /d.o.f ~ 1.15). 

Now we see that with this modification of the expression of the size of the front, the results 
for the cumulants shown in figure 4.12 (full lines) are in excellent agreement with the numerical 
data over the whole range of N. 
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Chapter 5 

Spatial correlations 



So far, all models and calculations aimed at describing QCD scattering amplitudes assumed uni- 
formity in impact-parameter space, or, decoupling of the evolution between different points in the 
transverse space. Indeed, we considered one- dimensional models while to fully describe the impact 
parameter, two dimensions are necessary. We shall address here the issue of the correlations of 
the QCD evolution between different impact parameters. 

Contents 

5.1 Relevance of one-dimensional models 

5.1.1 A model incorporating an impact-parameter dependence 

5.1.2 Numerical evaluation of the correlators 

5.2 Computing the correlations 

5.2.1 Basic features of the model 

5.2.2 Formulation of the calculation of the correlations .... 

5.2.3 Analytical expression for the correlations 

5.2.4 Numerical simulations 



EH 
EH 
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ED 
ED 
EH 



5.1 Relevance of one-dimensional models 

So far, we have argued that high-energy scattering in QCD at fixed coupling and fixed impact 
parameter is in the universality class of the stochastic FKPP equation (Chap. [2|, which is an 
equation with one evolution variable (time or rapidity in QCD), and one spatial dimension (x 
generically, or lnfc 2 ~ ln(l/r 2 ) in QCD). From the very beginning, we have simply discarded 
the impact parameter dependence. It is important to understand that the spatial variable and 
the impact parameter play different roles, and thus, the impact parameter may a priori not be 
accounted for by a two-dimensional extension of the FKPP equation. 

There are general arguments to support the assumption that the QCD evolution is local enough 
for the different impact parameters to decouple through the rapidity evolution, which we are now 
going to present. 

Let us start with a single dipole at rest, and bring it gradually to a higher rapidity. As was 
explained in Chap. [2] during this process, this dipole may be replaced by two new dipoles, which 
themselves may split, and so on, eventually producing a chain of dipoles. Figure [272] pictures one 
realization of such a chain. 



According to the splitting rate given in Eq. (2.1 1, splittings to smaller-size dipoles are favored 



and thus, one expects that the sizes of the dipoles get smaller on the average, and that in turn, the 
successive splittings become more local. The dipoles around region "1" and those around region 



"2" should have an independent evolution beyond the stage pictured in Fig. 2.2 Further splittings 
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will not mix in impact parameter space, and thus, the traveling waves around these regions should 
be uncorrelated. For a dipole in region 1 of size r to migrate to region 2, it should first split 
into a dipole whose size is of the order of the distance A6 between regions 1 and 2, up to some 
multiplicative factor of order 1. (We assume in this discussion that the dipoles in region 2 relevant 
to the propagation of the local traveling waves, that is, those which are in the bulk of the wave 
front, also have sizes of order r). Roughly speaking, the rate of such splittings may be estimated 



from the dipole splitting probability (2.1 1: it is of order a[r 2 j (A6) 2 ) 2 , while the rate of splittings 
of the same dipole into a dipole of similar size in region 1 is of order a. Thus the first process is 
strongly suppressed as soon as regions 1 and 2 are more distant than a few units of r. Note that 
for Ab > 1/Q S , saturation may further reduce the emission of the first, large, dipole leading to an 
even stronger suppression of the estimated rate. 

What could also happen is that some larger dipole has, by chance, one of its endpoints tuned 
to the vicinity of the coordinate one is looking at (at a distance which is at most |Ar| <C 1/Q s (y)), 
and easily produces a large number of dipoles there. In this case, the position of the traveling wave 
at that impact parameter would suddenly jump. If such events were frequent enough, then they 
would modify the average wave velocity and thus the one-dimensional sFKPP picture. We may 
give a rough estimate of the rate at which dipoles of size smaller than Ar are produced. Assuming 
local uniformity for the distribution n of the emitting dipoles, the rate (per unit of ay) of such 
events can be written 

' f1 ''" [ **«r*)(^ (5.1) 



r >Ar r J e<Ar " \ r o J 2ire 2 {r -e) 

where we integrate over large dipoles of size r > Ar emitting smaller dipoles (of size e < Ar) with 
a probability d 2 e r 2 . / (2ire 2 (ro — e) 2 ). The factor (e/ro) 2 accounts for the fact that one endpoint of 
the dipole of size ro has to be in a given region of size e in order to emit the dipoles at the right 
impact parameter. To estimate this expression, we first use n(ro) — T(ro)/a 2 and approximate T 

by 

T(r ) = 8(r - l/Q s ) + (r 2 Q 2 s )^ 8(1/Q S - r ). (5.2) 

The front is replaced by 1 above the saturation scale (for ro > 1/Q S ) and by an exponentially 
decaying tail for r < 1/Q S . Using r — s ~ r in the emission kernel, the integration is then easily 
performed and one finds a rate whose dominant term is 

tt (( Ar) 2 Q 2 )^ (53) 



2a 2 1 - 7 c 

For (Ar) 2 -C (a 2 ) 1 / 7c /Q 2 , i.e. ahead of the bulk of the front, this term is parametrically less than 
1 and is in fact of the order of the probability to find an object in this region that contributes to 
the normal evolution of the front |114| . Hence there is no extra contribution due to the fact that 
there are many dipoles around at different impact parameters. 

The arguments given here are based on estimates of average numbers of dipoles, on typical 
configurations, and we are not able to account analytically for the possible fluctuations. As we have 
seen through this review, the latter often play an important role. As a matter of fact, in the physics 
of disordered systems, rare events sometimes dominate. So before studying the phcnomcnological 
consequences of the statistical picture of high-energy QCD based on a one-dimensional equation, 
one should check more precisely the locality of the evolution in impact parameter. 

A numerical check was achieved in the case of a toy model that has an impact-parameter 
dependence in Ref. |117 . Let us briefly describe the model, before presenting the main numerical 
results. 

5.1.1 A model incorporating an impact-parameter dependence 

In order to arrive at a model that is tractable numerically, we only keep one transverse dimension 
instead of two in 3+1-dimensional QCD. However, we cannot consider genuine 2+1-dimensional 
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QCD because we do not wish to give up the logarithmic collinear singularities at X2 = xo and 
x i — x i- Moreover, QCD with one dimension less has very different properties at high energies 



118 . Starting from Eq. (2.1 1, a splitting rate which complies with our requirements is: 

dP 1 |x i| , . 

~m^ = ^~\ ji i dx 2- (5.4) 

d(ay) 4 |x 02 ||a:i2| 

We can further simplify this probability distribution by keeping only its collinear and infrared 
asymptotics (as in Ref. [119] ). If | £02! <C |xoi| ( or the symmetrical case \xi2\ <C | a^oi I) ■ the 



probability reduces to dx2/\xo2\ (^2/1^12! resp.). The result of the splitting is a small dipolc 
(xo,X2) together with one close in size to the parent. So for simplicity we will just add the small 
dipole to the system and leave the parent unchanged. In the infrared region, a dipole of size 
1^02 1 3> I a^oi I is emitted with a rate given by the large-|.To2| limit of the above probability. The 



probability laws (2.1 1,(5.4) imply that a second dipole of similar size should be produced while the 
parent dipole disappears. To retain a behavior as close as possible to that in the collinear limit, 
we will instead just generate a single large dipole and keep the parent. To do this consistently one 
must include a factor of two in the infrared splitting rate, so as not to modify the average rate of 
production of large dipoles. 

Let us focus first on the distribution of the sizes of the participating dipoles. (The simplifying 
assumptions made above enable one to choose the sizes and the impact parameters of the dipoles 
successively). We call r the modulus of the emitted dipole, ro the modulus of its parent and we 



define Y — ay. The splitting rate (5.4 1 reads, in this simplified model 



, v =0(r- r )— J- + 0(r - r) — , (5.5) 
dY r z r 

and the original parent dipole is kept. Logarithmic variables are the relevant ones here, so we 
introduce 

p = log 2 (l/r) or r = 2- p . (5.6) 
We can thus rewrite the dipole creation rate as 
dP ^ 

= e(p - p) 2"-p« log 2 dp + 6{p - p ) log 2 dp. (5.7) 

To further simplify the model, we discretise the dipole sizes in such a way that p is now an integer. 
This amounts to restricting the dipole sizes to negative integer powers of 2. The probability that 
a dipole at lattice site i (i.e. a dipole of size 2~ l ) creates a new dipole at lattice site j is 

dPi^j _ r +1 dP p ^ p _ flog 2 j > i 
dY -J pj dY j<i (5 ' 8) 

The rates dPi± /dY for a dipole at lattice site i to split to any lattice site j > i or j < i respectively 
are then given by 

dPi+ dPi-jfj . .. dPi- dP^j ■ 

j=» 3=0 

where we have restricted the lattice to < i < L, for obvious reasons related to the numerical 
implementation. 

Now we have to address the question of the impact parameter of the emitted dipole. In QCD, 
the collinear dipoles are produced near the endpoints of the parent dipoles. Let us take a parent 
of size ro at impact parameter 60 . We set the emitted dipole (size r) at the impact parameter b 
such that 

b = b ± r -^^ (5.10) 
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where s has uniform probability between and 1. It is introduced to obtain a continuous distri- 
bution of the impact parameter unaffected by the discretisation of r. This prescription is quite 
arbitrary in its details, but the latter do not influence significantly the physical observables. Each 
of the two signs that appear in the above expression is chosen to be either + or — with equal 
weights. We apply the same prescription when the emitted dipole is larger than its parent. 



Now that we have introduced a branching process similar to QCD dipole evolution, we must 



define the scattering amplitude. We have explained above (see Sec. 2.1 ) that in QCD, the scattering 
amplitude of an elementary probe dipole of size r t — with a dipole in an evolved Fock state 
is proportional to the number of objects which have a size of the same order of magnitude and 
which sit in a region of size of order around the impact point of the probe dipole. Since in our 
case, the sizes are discrete, the amplitude is just given, up to a factor, by the number of dipoles 
that are exactly in the same bin of size as the probe, namely 

T(i, bo) = a 2 s x #{dipoles of size 2~ J at impact parameter b satisfying | fo — £> 1 < r i/2}- (5-11) 

Finally, we must enforce unitarity, that is, the condition 

T(i,b)<l (5.12) 

for any i and b. This condition is expected to hold due to gluon saturation in QCD. However, 
saturation is not included in the original dipole model. The simplest choice is to veto splittings 
that would locally drive the amplitude to values larger than 1. In practice, for each splitting that 
gives birth to a new dipole of size i at impact parameter b, we compute T(i, b) and T(i, b ± r^/2), 
and throw away the produced dipole whenever one of these numbers gets larger than one. 

Given the definition of the amplitude T, this saturation rule implies that there is a maximum 
number of objects in each bin of size and at each impact parameter, which is equal to -/V sat = l/a 2 . 



5.1.2 Numerical evaluation of the correlators 

We have implemented this model numerically. Let us discuss how we operated this implementation 
and the results we obtained. 

We take as an initial condition a number N sat of dipoles of size 1 (i = Q) , uniformly distributed 
in impact parameter between — ro/2 and ro/2. The impact parameters bj that are considered are 
respectively 0, 10 -6 , 10~ 4 , 10~ 2 and 10 _1 . The number of events generated is typically 10 4 ,which 
allows one to measure the mean and variance of the position of the traveling waves to a sufficient 
accuracy. 

We have checked that at each impact parameter, we get traveling waves whose positions grow 
linearly with rapidity at a velocity less than the expected mean-field velocity for this model. iV sa t 
was varied from 10 to 200. 

Figure |5.1| represents the correlations between the positions of the wave fronts at different 
impact parameters, defined as 

(p s (Y,h)p s (Y,b 2 )) ~ (pB<y,h))<P.(y,l>2))- (5-13) 

We set -/V sat to 25 in that figure, but we also repeated the analysis for different values of iV sa t 
between 10 and 200. 

We see very clearly the successive decouplings of the different impact parameters in Fig. |5.1[ 
from the most distant to the closest one, as rapidity increases. Indeed, the correlation functions 
flatten after some given rapidity depending on the difference in the probed impact parameters, 
which means that the evolutions decouple. This decoupling is expected as soon as the traveling 
wave front reaches dipole sizes which are smaller than the distance between the probed impact 
parameters, i.e. at Y such that 1 1>2 — &i| ~ l/Qs(Y) = 2~ Ps ( Y \ (We shall further comment on 
this decoupling in the next section). From the data for p s (Y), we can estimate quantitatively 
the values of the rapidities at which the traveling waves decouple between the different impact 
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{ Ps (Y, 0) Ps (Y, b)) - ( Ps (Y, 0))( PS (Y, b)) 



14 




5 



Y 



Figure 5.1: Correlations of the positions of the traveling wave fronts between different impact 
parameters in the toy model of Sec. |5.1| The points where the correlations flatten correspond to 
the decoupling of the waves in the corresponding regions of impact parameter. 



parameters. (It is enough to invert the above formula for the relevant values of 62 — bi). These 
rapidities are denoted by a cross in Fig. |5. 1| for the considered impact parameter differences. Our 
numerical results for the correlations are nicely consistent with this estimate, since the correlations 
start to saturate to a constant value precisely on the right of each such cross. 

We conclude that the different impact parameters indeed decouple, as was expected from a 
naive analytical estimate. What is true for our toy model should go over to full QCD, since we have 
included the main features of QCD. When looking at the data more carefully however, it turns 
out that the model with impact parameter does not reduce exactly to a supposedly equivalent 
one-dimensional model of the sFKPP type. This is a point that would deserve more work. We 
refer the reader to Ref. 117 for all details of our numerical investigations. 

However, even if one takes the statistical decoupling of impact parameters as soon as Ab > 
l/Q s (y) for granted, there may still be some effective correlations persisting at large rapidities 
since two points in impact-parameter space share some common history. Indeed, fluctuations need 
some rapidity to affect saturation scales, and history may be remembered way after the rapidity at 
which the decoupling happened. Such effects are negligible at small rapidities, as was just shown 
numerically, but may be crucial at large rapidities. We shall now investigate this point. 



5.2 Computing the correlations 

Our method will consist in proposing a simple toy model which contains the main physical features 
of QCD, which may be implemented as a Monte-Carlo event generator and for which analytical 
calculations will be possible. In these respects, our approach follows the one developed in the 
previous section, but while the latter work was purely numerical, our main results will consist 
in analytical expressions of the correlation of the saturation scale between two points in impact- 
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parameter space, as a function of the distance between the points and as a function of the rapidity. 



The content of this section was published in Rcf. 120 121 . 

The model that we will study will have the following characteristics. With respect to QCD, 
we assume the following simplifications: (i) Dipoles evolve by giving birth to one dipole of half 
size (the left or the right half of the parent dipole) , or to one dipole of double size (in such a way 
that the parent be the left or right half of its offspring) at some fixed rates, (ii) dipoles do not 
disappear in the evolution, that is to say, the parent dipoles are not removed, (Hi) the positions 
and dipole sizes are discrete, and (iv) the configuration space of the dipoles is a line instead of the 
full two-dimensional space. We thus give up two main properties of the QCD dipole model: The 
collinear singularities, which cause the dipole endpoints to emit an arbitrary number of dipoles 
of arbitrarily small sizes, and the continuous and two-dimensional nature of the dipole sizes and 
positions. The first simplification is the diffusion approximation, which has been studied in the 
context of BFKL physics (see e.g. Ref. |l22|), but which was not assumed in the previous model. 
The second simplification was instead already assumed in there. These model simplifications may 
introduce some artefacts, but that we believe are under control, and many results which we will 
obtain within such simple models are likely to apply to QCD since they will not depend on the 
details. 

Let us now specify completely the model. According to the evolution rules given above, starting 
from a dipole of size 1, the sizes of all dipoles present in the system after evolution are powers of 
2. In practice, we shall only consider fractions of 1, i.e. the sizes may be written as 2 1_fc , where 
k > 1. For each value of fc, there are 2 k ~ 1 possible values of the position b of the center of the 
dipoles: b = —\ + 2~ k , — \ + 3 x 2~ fe , • • • , | — 3 x 2~ k , \ — 2~ k . Let us number these bins by the 
index < j < 2 k - 1 - 1 running from the negative to the positive positions. The model may be 
represented as a hierarchy of bins that contain a discrete number of dipoles, see Fig.|5.2| Note that 



to any given impact parameter b between — i and | corresponds one unique bin at each level of size. 
For example, at position b — — |, one sees the bins (fc = l,j = 0), (k = 2,j = 0), (k = 3, j = 0) 
etc... At position —0.2, one sees the bins (fc = l, j = 0), (k = 2, j = 0), (fc = 3,j = 1) etc... More 
generally, at position — \+ A6, one sees (fc, [Ab x 2 k - 1 ]), where the square brackets represent the 
integer part. 

During the rapidity (or time) interval dt, a dipole in the bin (fc, j) has a probability adt to give 
birth to a dipole in the bin (fc + 1, 2j), adt to give a dipole in the bin (fc + 1, 2j + 1), and f3dt/2 
to give a dipole in the bin (fc — l,j/2) if j is even and (fc — 1, (j — l)/2) if j is odd. Note that dt 
may be infinitesimal (which is generally speaking convenient for analytical calculations), but also 
finite (which is convenient for numerical simulations). 

As for the saturation mechanism, we assume the simplest one: We veto splittings to bins which 
already host the number N of dipoles. 

We can consider that the number density of "gluons" of a given size seen at one impact pa- 
rameter is proportional to the number of dipoles in the corresponding bin (fc, j). As rapidity is 
increased, the occupation of the bins with low values of fc gets higher until the number of objects 
they contain reaches TV. The subsequent filling of the bins indexed by larger values of fc (smaller 
dipole sizes) can be seen as the propagation of traveling wave fronts at each impact parameter, 
with possibly complicated relationships between them. The (logarithm of the) saturation scale 
X(b,t) at impact parameter b is related to the position of the front seen there at time t. There 
are several equivalent ways to define the position of the front. It could be, for example, the largest 
value of fc for which the number of objects becomes some given fraction of N . (Later, we will use 
a slightly different definition) . 



5.2.1 Basic features of the model 

Let us denote by ri(k,j){t) the number of dipoles present in the bin (fc, j) at time t. Then, according 
to the rules given above, we can write the following stochastic evolution equation: 



n (k,j){t + dt) 



N,n 



(*) 



<*(Jfc-l,[j/2])(*) 



5 P/2 



(t) + 5, 



13/2 

0+1,2.7+1) 



«(*) 



(5.14) 
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Figure 5.2: The hierarchical structure of the model. Each box represents a bin which may contain 
up to N dipoles of given sizes (vertical axis) and positions in impact-parameter space (horizontal 
axis). The conventional numbering of the bins that we have chosen is also shown for k = 1, 2, 3. 



where the 5f k ^ are drawn according to the binomial distribution 



Proba 



(5.15) 



This is a rather complicated equation which we do not know how to solve except numerically. 

This model does not a priori look like a stochastic FKPP model. We may assume uniformity 
in impact parameter: This would amount to imposing the same S a and 5^/ 2 respectively for all j 
at any given k. In this case, the model would be projected to the FKPP class, but by definition, 
this would wash out the fluctuations between the different impact parameters. This simplified 
model, that we call "FIP" (for "Fixed Impact Parameter", since effectively, the model is completely 
defined by a single impact parameter) is nevertheless useful since it provides a benchmark to 
evaluate how the fluctuations between different impact parameters may alter the FKPP picture. 
In this paper, we will rely on (and check again in the case of our model) the conclusion reached in 
Ref. 1117 that thanks to saturation, locally at each impact parameter, the full model is still well- 



described by a one-dimensional FKPP equation, and the fluctuations between different positions 
in impact-parameter space do not qualitatively change the picture. 

Let us first apply the treatment of FKPP equations exposed in Chap. [4] to the FIP case. We 
know that the large-rapidity realizations of the model are stochastic traveling waves, whose main 
features can be determined from a simple analysis of the linear part of the evolution equation. In 
this model, only the number of dipoles n k in the bins say (fc,0) (i.e. at impact parameter — |) is 
relevant. The evolution equation readfQ 



rife it + dt) 



N,n k (t) + 5^(1:) 



Ci(*) 



(5.16) 



The mean-field (or Balitsky-Kovchegov) approximation to the evolution leads to the equation 



n k (t + dt) = min 



N, n k (t) + adtnk-i(t) + f3dtn k+1 (t) 



(5.17) 



could also write 25"^'^ (t) instead of the last term in Eq. ( 5.16 I. (This may even be a more literal implemen- 
tation of the FIP approximation). But this would not make a large difference, which anyway, we would be unable 
to capture analytically. 
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where the n k are now real functions of k. The linearized equation (equivalent to the BFKL 
equation) is simply obtained by discarding the "min" in the previous equation: 

n k (t + dt) = n k (t) + adtn k _ 1 (t) + f3dtn k+1 (t). (5.18) 

From standard arguments, we know that for asymptotically large t and TV, the velocity of the wave 
front, that is the time derivative of the position X(t) of the front, is given by 

dX f, » . . 

v ° = dF = (7c) ' ( } 



where w(7) is the eigenvalue of the kernel of the linearized evolution equation (5.18) corresponding 
to the eigenfunction e _7fc , namely 

= -r In (1 + adt e 7 + f3dt e~ 7 ) , (5.20) 

dt 



and 7 C minimizes lu( , -/)/j. We recall that dt may be finite or infinitesimal, in which case Eq. (5.201 
is to be understood as the derivative of ln(l + •) at the origin. 

Our aim is to study the corre latio ns between the point at position b = — \ in transverse space 
(left edge of the system, see Fig. 5.2 1 and the one at position b = — | + Ab with < Ab < 1. We 



calculate the average of the squared difference of the positions of the front between these points, 
which is formally related to the two-point correlation function of (the logarithm of) the saturation 
scales, and which we deem a good estimator of the spatial fluctuations of the saturation scale. 
In the hierarchical model, all bins with index k less than or equal to k^ = 1 + [— log 2 Ab] (the 
notation "[• • •]" stands for the integer part) and j ' = overlap both impact parameters, and thus 
the dipoles of size larger than 2~ kAh seen at these points are exactly the same. For k > k^b 
instead, the bins seen at the two points are distinct and nonoverlaping. So in particular, in our 
model with j3 = 0, as soon as the position of the front at one point or at the other is larger than 
/cAb, that is to say, as soon as there are of the order of N dipoles in the bin (k&b,j = 0), then the 
evolutions are completely uncorrelated at the two points in the corresponding bins. (We expect 
that for finite j3 of order 1, the discussion would not be qualitatively changed.) This matches to 
the picture that we may infer for the QCD dipole model: The dipoles at two positions in impact- 
parameter space separated by a distance larger than the typical saturation scales in that region 
evolve (almost) independently towards larger rapidities. Note that choosing pairs of points around 
impact parameter 0, one with positive impact parameter and another one with negative impact 
parameter, would not satisfy this property, due to the rigidity of the sizes and positions of the 
dipoles. Indeed, these two points would decorrelate very soon in the evolution since their common 



ancestors necessarily sit in the bin (k = 1, j = 0), see Fig. 5.2 

As a consequence of these features of QCD reproduced in the toy model, studying two-point 
correlations between points in impact-parameter space as a function of their distance Ab and of 
the time (=rapidity) t is equivalent to studying the time dependence of the correlations of the 
saturation scales of two realizations of the model whose evolutions are identical until the tip of the 
front reaches k/±b- On the average, it takes a time tA& = (&A6 — v being the mean velocity 

of the individual fronts, for the front whose tip is at k&b = 1 at the beginning of the evolution 
to have its tip at /ca&- Then the bins such that k > k&b evolve independently between the two 
realizations over the remaining time interval 

At = t- t Ab , with t Ab = — 1<3g2 A ^ . (5.21) 

v 

Note that this is very close to assuming that the realizations are identical for t < t&b and completely 
uncorrelated for t > t^b- 

From this discussion, we see that the basic input of our calculation will be the mechanism for 
the propagation of a FKPP front which was explained in Chap. |4j We will review it in the next 
subsection, then we will proceed to the formulation of the calculation of the correlations. 
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5.2.2 Formulation of the calculation of the correlations 



In line with the above discussion, we wish to compute the correlations of the position of two fronts 
whose evolutions are identical for t < t/\b and uncorrelated for t > £a6- Note that strictly speaking, 
we would need to keep the content of all bins k < k^b identical between the two realizations at all 
times, even after time t&b- But these two formulations give quantitatively similar results. 

Let us introduce X(t , t) the position of the front at time t in the frame in which A(t , to) = 0. 
We focus on what happens slightly before the initial time t a . According to the mechanism of 
front propagation explained in Chap. [4] on one hand, X(to — dto,t) = X(to,t) + Ubd^O if no 
fluctuation has occurred between times to — dto and to, on the other hand, X{to — dto,t) = 
X(to,t) + wbd^o + R{t — to, 5) if a fluctuation has occurred at a position S ahead of the front 



(which happens with probability p(S)dS dto). «bd was defined in Eq. (4.521, R(t,S) in Eq. (4.571 
and p(S) in Eq. ( 4.54[ ). It is straightforward to write an equation for the generating function of 
the cumulants of X: 



d 

dt 



i n ( e A*(*o,t; 



At>BD + dS 



>p(<J) (< 



\R(t-t ,S) 



1 



(5.22) 



One now considers two such independent fronts and add up the generating functions. One gets 

(5.23) 



d 
di 



- In ((eAXxCto,^ ^ e -Ax a (to,t)^ = / d5p{5) ( 



e XR(t-t a ,S) _|_ e -XR(t-t ,S) _ 2 



Expanding for A close to 0, the coefficients of the second power of A obey the equation 

d 



dt 



(pd - X 2 f) = 2 / d6p{6)R 2 (t - t Q , S), 



(5.24) 



where we have used the fact that X\ and X 2 are independent random variables for t > to, and we 
have traded t for t in the derivative, taking advantage of the fact that both Xi — X 2 and R only 
depend on t — t . In practice, t will be equal to £a6j the time at which the tip of the single front 



reaches k^b- From Eq. (5.211, this time is [— log 2 Ab]/v 



We see that the basic ingredient is the time evolution of the shift of the front due to a forward 
fluctuation. This shift was given in Eq. (4.571. It involved the expression for ipg m Eq. (4.461. 



In order to write a compact expression for R, it is interesting to note that ips is related to some 
Jacobi d function 1123 . Since 



d A {z\q) = l + 2^(-l)"cos(2nz)g™ 2 , 



(5.25) 



we may rewrite Eq. (4.461 as 



The notation 



n(a + p) 



q -0 4 



?r(a - p) 



= P. 4 = 



q = e 



= e 



(5.26) 



(5.27) 



has been introduced. Using Eq. ( |4.57 1 and performing the appropriate expansion for small a and 
a, we arrive at an expression for R(t, S) in terms of the 04-function which is particularly compact: 



R(t, S) = — In 

7c 



with 



+ CO 



^ 4 (0k) = 2^(-l)" +1 nV _1 . 



(5.28) 



(5.29) 
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It is actually quite natural that the Jacobi theta functions appear, since the latter are defined as 
solutions of the one-dimensional heat equation with periodic boundary conditions. 

We turn to the analysis of the obtained result. First, for large y, only the fundamental mode 
contributes significantly to ipg. Looking back at Eq. (4.461, we see that higher harmonics would 
give a series of exponentially decreasing corrections. But at a finite time, a large number of modes 
have to be taken into account, typically all modes such that n < (L / n) yj2 / lu" ("f c )t. A few low- 
lying modes are not enough to describe the small-time behavior. Instead, it is a saddle point (in 
an appropriate integral reformulation) that dominates the sum ( 4.46| . In this regime, it would 
be useful to find a way to write the series of harmonics such that at asymptotically large y, only 
the first term contributes instead of the whole series. This is actually possible using the Poisson 
summation formula 

+00 +00 

E /(»)= E / dxf(x)e- 2i * k *. (5.30) 



n— — 00 



k=-c 



In order to get R(t, 5), we need the value of ip$ at p — a. Hence we choose 

e 7c<5 2 
f(x) = — sin nxa sin nxa q x ~ 1 e t7TX . 



(5.31) 



We then perform the integral over x in the r.h.s. of Eq. (5.301. Introducing 7+ = a + a and 
7_ = a — a, we get the following expression for ips- 



tps(y,a) 



alaS I 



4L 2 q y — In q 



+00 

E 

k=— 00 



(2fc-l+-y+r 

e 41n " 



(2fc-i-7+r 

e 41n a 



(2k-l + -i_y 



(2fc-l— 

e 41n ' 



(5.32) 

Since we eventually want to apply Eq. (4.571 in order to get an expression of the shift of the front, 
we expand the latter formula for a, a <C 1. The leading order reads 

a~fc$ 



T q(-\nqY/ 2 ~L? 



2- +°° 
7r aa 



[n 2 (2k- l) 2 + 2\nq] 



e 41n 9 



(5.33) 



k=l 



The shift of the front due to a fluctuation is obtained from ips with the help of Eq. (4.571: 

+00 



R(t, 6) = — In i 1 + C 2 

7c 



/2L 2 e^ 5 

(W( 7c )t)5/2 



E 

k=l 



(2k -iy - 



2 W"(7c)* 



L 2 



g 2u"( 1c )l 



(5.34) 



q is the function of t given by Eq. (5.27 1. This formula is extremely useful, since the series indexed 
by k converges fast. Even for moderately large values of t, a few terms accurately describe the 
whole function. This is actually the best formula for numerical evaluations of R. 

We shall now examine the limit of small t (y -c 1). Then only the term k = 1 has to be kept. 
The expression for R boils down to 



1 



R(t, 5) = — In 1 + C 2 



y/2e / " 5 L 2 



'(7TU/'( 7C )) 5 / 2 i 5 /2 



(5.35) 



As a final remark, let us note that the Poisson summation (5.30 1 that we have used to rewrite 
the series of harmonics corresponds to a Jacobi identity for the i9 functions 123 . Equation (5.341 
results from Eq. (5.281 with the replacement 

^ 1 



+ OO 



0,A*(O|g) 



2 q(-\nqfl 2 



( 7I " 2 (2fc - I) 2 + 2 Inq) e «W 



(5.36) 



fc=i 



72 



CHAPTER 5. SPATIAL CORRELATIONS 



5.2.3 Analytical expression for the correlations 

With the elements presented in the previous sections, we can write the expression for of 2 = 
((Xi — X2) 2 ). It is enough to insert the expression for the probability of fluctuations (Eq. (4.54)) 
and for the time-dependent shift (Eq. (5.28 I ) into Eq. (5.241: 



dcr\ 2 
dt 



2d 



+00 



dS e" 1 ^ In 2 



(5.37) 



where for 9 9 i?4(0|g) we use either one of the equivalent expressions (5.29), (5.361 according to the 
limit that we want to investigate. We now have to fix the value of L. In Rcf. 1 1 14 , L was taken 
to be a constant. (The phenomenological model predicted L = L$ = lniV/ 7c , but empirically, we 
saw that it was better to add a subdominant correction, namely L = Lq + — lnL + const.) In 
this case, a change of variable can be made in the integrand. All the parameters may be factored 
out, leaving us with a simple numerical integral to perform: 



+00 



dx 



ln 2 (l + a;) = 2C(2) 



(5.38) 



Thus 



dt 



TT 2 dC 2 



-d q M0\q)} 



(5.39) 



3 7 3 L 3 

Replacing the product of the unknown constants by Eq. ( |4.65| ) and q by Eq. (5.271 and integrating 
over the time variable between and At = t — t^, we arrive at a parameter-free expression for 
a 2 2 as a function of At, namely 



2 

a 12 



2tt 2 
3 7 3 L 



dq 

q 



-W<%)] 



(5.40) 



We now investigate the two interesting limits, i.e. At 3> L 2 and At C L 2 . For large At, the 
integral is dominated by the region q — > 0, thus — c? 9 t?4(0|g) may be replaced by its value at q = 
(— <9 g $4(0|0) = 2). Performing the remaining integration, we get 



At»L 2 



27TV'( 7c ) 

3 7 3 L 3 



At, 



(5.41) 



which is twice the second-order cumulant of the position of the front in a fixed impact-parameter 
model, see Eq. ( |4.66 ). For small At instead, say L <C At <C L 2 , we use the expansion of dqtii^q) 
for q —¥ 1, i.e. the first term in Eq. (5.361, which reads 



d q M0\q) = - 



Vvr 7T 2 + 2 In q 
~T q{-lnq) 5 /< " 



g41n, 



(5.42) 



Equation (5.40) boils down to the following expression: 



2 

CT 12 



At<SL 2 




2tt 3 



w"( 7c )Ai 



exp 



L 2 



2o;"( 7c )At 



(5.43) 



So far, we have chosen the size of the front L constant, of the order of Lq. Another possible 



model for L would be to promote it to a function of 5 at the level of Eq. (5.371, namely 



L = Lq + S + const, 



(5.44) 



where the constant has to be determined empirically. This choice takes maybe into account more 
accurately the extension of the front by S generated by the fluctuations. The (^-integral cannot 



be performed analytically in Eq. (5.371 except in some limits, so a priori, there is no simpler 
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expression than Eq. ( 5.37| . Thus we need to know the values of C\ and C2 individually. We can 
consider that C\ = 7 C is the natural normalization of the probability distribution p(5). Then, we 
must set C2 = tt 2 lj" (7c) /7c i n order to satisfy Eq. (4.651. 

The above-mentioned two models, in which L is either constant or <5-dependent, differ by 
subleading terms in the large- L limit. Since the values of 5 which dominate the ^-integral in 
Eq. (5.371 are of order — InLo, like the first correction to Lq in the case of constant L, the models 



are not expected to differ significantly. We will check this statement numerically. 



Scaling 



Looking back at Eq. (5.401, we see that a\ 2 has a nice scaling property. Indeed, we may rewrite 
the latter equation as 



D 



'12 



7c(w - v) y e -ToC«o-«)A< q 



(5.45) 



in terms of the properties of a single front (its velocity v and the diffusion constant D whose 
analytical expressions were given in Eq. (4.521 and (4.66 1), where Vo can be read in Eq. (5.191. In 
particular, we have the following scaling: 



"12 

DAt 



= function[(w — v)At]. 



(5.46) 



From Eq. (5.431, we see that the function in the right-hand side is exponentially damped when its 
argument is smaller than 1, i.e. parametrically for At <C L 2 . 

Once one knows the characteristics of the traveling waves in the FIP model (i.e. v and D), 
this scaling of the correlations is a pure prediction. Thus it will be interesting to check it in the 
numerical calculations. 



Limits on the validity of the calculations 

Let us try and evaluate the limits on the validity of our calculations. The latter were essentially 
based on the assumption that the eigenvalue 7 = 7 C of the kernel w dominates. While this 
statement is clearly true at large times, when the traveling-wave front is well formed, it must 
break down at early times right after a fluctuation has occurred: Indeed, a fluctuation has an 
initial shape that is far from the one of the asymptotic front, see Fig. 4.10 

We wish to estimate the order of magnitude of the dispersion of the relevant eigenvalues about 
7c. To this aim, neglecting for the moment the boundary conditions and the prefactors, we write 
the solution of the general branching diffusion as 



u(At,k) 



drf e 



-7fc+td(7)At 



(5.47) 



The interesting values of k are the ones around the position of the wave front, therefore we write 
k = voAt + 5k, where 5k is of the order of the size L of the front. Expanding u>(-f) about j c , we 
write 



u(At, v At + 5k) ~ e 



(5. 



where 5^/ = 7 — 7 C . It is clear from this equation that the relevant values of £7 are of the order 
of 5k/ (w"(7c)Ai). Since the order of magnitude of 5k is the size L of the front, we would a priori 
conclude that the dispersion of 7 around 7 C is small and hence that the calculation is valid as soon 
as At > L. 

However, we have also expanded 07(7) to second order. This means that for a generic kernel 
w, we have neglected terms of the form ^oj^(j c )(5j) 3 At ~ L 3 /(At) 2 (which would fit in the dots 



in Eq. ( 5.48 1) . The expansion is a good approximation if the latter term is small, i.e. if 

At > L 3/2 . 



(5.49) 
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Back to impact-parameter space 

So far, we have been working with the minimal model, consisting in two realizations of the FIP 
model which evolve in the same way until their common tip reaches k^, and which decorrelate 
for k > fcA6- The only relevant parameter which determined the decorrelation of the positions of 
the fronts of the realizations was the time At = t — t&b after the tip had reached k^. We now 
wish to discuss the transcription of the obtained results to impact-parameter space, which was 
our initial problem. 



To this aim, we will of course make use of Eq. (5.21 I to express with the help of the mean 
front velocity v. But we also need a length scale to which the distance in impact-parameter space 
Ab may be compared. The natural length is the dipole size at the position of the front, namely 



l.(t) 



-X(t) 



l s {t 



Ab 



-vAt 



On the other hand, according to Eq. (5.21 I and disregarding the integer part operator 
kAb and the tip of the front is ahead of the bulk X(t&b) by L: k&b 
previous equation, we may now express At as a function of Ab and of the length scale l s (t): 



(5.50) 

log 2 Ab = 
X(tAb) + L. Using the 



At = - 
v 



L + log 2 



Ab 



l s (t) 



The scaling (5.461 reads 



L + log 



Ah 

2 hit) 



'12 



I? 



function 



L + log 



Ab 

2 h{t) 



L 2 



(5.51) 



(5.52) 



This formula, together with the behavior of the scaling function (see Eq. ( |5.43[ )), shows that there 
is little 6-dependence until \og 2 (Ab/l s (t)) ~ L 2 , that is to say, until Ab ~ l s (t)e' 



const X L 



In other 

terms, the size Ab of the domain around impact parameter b in which the fluctuations in the 
position of the fronts are negligible is, in notations more familiar to QCD experts, 



Ab' 



O constxln 2 (l/o:^) 

QAb) : 



(5.53) 



where Q s {b) is the usual saturation momentum at impact parameter b. Note that since the fronts 
are statistically independent as soon as Ab x Q s (b) > 1, this result may seem a bit surprising: 
It says that the effective correlation length between different points in impact-parameter space is 
much larger than 1/Q s {b) in the parametrical limit of small a s . 



5.2.4 Numerical simulations 

In this section, we confront our analytical calculations to numerical simulations of the toy model. 
First, we consider the full model and test the validity of the assumption that the minimal model is 
a good approximation to the full model also for j3 ~ 1, i.e. when splittings to larger-size dipoles are 
authorized. Second, we compare the minimal model to the analytical results for the fluctuations 



between different positions in impact-parameter space (given essentially by Eqs. (5.371,(5.40)) 



Full model 



The model defined by Eqs. (5.14) and (5.151 is straightforward to implement numerically in the 



form of a Monte-Carlo event generator. The simplest is to store the number of dipoles in each 
bin in an array whose index i is related to k and j through i = 2 fc_1 + j. The splitting dynamics 
relates bin i to 2i (down left), 2i + 1 (down right) and [i/2] (up; the square brackets stand once 
again for the integer part). 

We have to deal with an array whose size grows exponentially with time. It is thus very difficult 
to pick large values of t, and thus also large values of N. Indeed, the relevant time scale grows 
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Figure 5.3: One event of the full model with a = (5 = 1, N = 100 and t — iA& = 4. Only the 
bins k > 5 are represented. (The bins for k < 5 all contain N dipoles.) The number of dipoles in 
each bin is proportional to the blackness which is displayed. We see that in the transition region 
close to blackness, nearby bins are often of similar grey levels, which illustrates the statement that 
the density of gluons varies significantly only over scales which are larger than the relevant length 
scale l s (t) (see Eq. ( |5.50| |). 



with TV like In 2 N, and consequently the minimum number of entries in the array one wants to 
consider grows like e N . In practice, we limit ourselves to t < 4 and N < 100. As for the time 
step dt, the most convenient is to take it small but finite. We set dt = 10 -2 . 

We start with one particle and evolve it for a few hundred units of time using the FIP version 
of the model. We obtain a traveling wave front, whose tip we eventually label k = 1. (The 
complete front sits in the bins k < 1). From the initial condition built in this way, we evolve all 
bins for which k < 1 using the FIP model, and all bins for k > 1 using the full model. One event 
is shown in Fig. |5.3| Although N and t are small in this calculation, we see that the regions in 
impact-parameter space which have similar numbers of dipoles are larger than the local length 
scale l s (t) (see Eq. ( |5.50| )). 



After the evolution times t = 3 and t = 4 respectively, we measure the position of the front at 
various impact parameters on a uniform tight grid ranging from — i to +|. We use the following 
definition of the position of the front: 

X(Ab,t) = k + 2^ — — — ^ — (5.54) 

k=k + l 

where ko is the largest k for which ^(/ c .[A&x2 fc - 1 ]) = N. Note that in principle, we could have 
chosen X(Ab,t) — ko. In practice however, because of the discreteness of k in our model, this 
choice would introduce artefacts which we do not expect in real QCD. 

We compute the squared difference of the front positions between the impact parameters — | 
and — | + Ab, and average over events. We plot the result as a function of t + log 2 Ab/v, where v 
is the average front velocity measured at impact parameter — | . 

We compare the results to the correlations obtained in the minimal model, i.e. when we consider 
two independent realizations of an initial front. We do not attempt to compare to our analytical 
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N = 100, P = 
0.6 , , , , , r 




0.5 1 1.5 2 2.5 3 3.5 4 



t — t&b or t + log 2 Ab/v 



Figure 5.4: <j\ 2 = ((Aj — A 2 ) 2 ) as a function of At = t — t^ in the full model with a = 1 and 
P = (lines with steps; one corresponds to an evolution time t = 3, the other one to t = 4) and 
in the minimal model. In the full model, £a/> = [— log 2 Ab]/v, where v is the measured velocity 
of the front at impact parameter — |. In the FIP model, t^b is a fixed time, and corresponds to 
the time at which the tip of the front reaches k^, the bin after which two uncorrelated evolutions 
take place. 



formulas since the values of N that we are able to reach are too small for the approximations that 
we had to assume to be relevant. 

The corresponding plot is displayed in Fig. |5.4| for N = 100, a = 1, (3 = 0, and in Fig. |5.5| with 
the same parameters except (3 = 2. First, we see that in the full model, the graph of a\ 2 exhibits 
steps, i.e. cr 2 2 is constant by parts. This is related to the hierarchical structure of the model: The 
correlations between b = — \ and any of the points at b > are identical; The same is true for 
—0.25 < b < 0, —0.375 < b < —0.25 etc... The logarithmic &-scale on the i-axis makes the widths 
of the steps all equal. Next, we see that for small t — t&b (i-e. impact parameters close to — |) 
there are very little fluctuations in the front positions. 

Finally, we see that for /3 = 0, as anticipated, the full model and the minimal ones coincide 
almost perfectly (Fig. |5.4| . For j3 = 2, i.e. when splittings towards larger dipole sizes are switched 
on and therefore new correlations appear beyond the ones taken into account in the minimal model, 
there are some quantitative differences for large t (Fig. 5.5). But we see that using the minimal 
model instead of the full model that keeps all impact parameters is a good approximation. This 
corroborates the conclusions of the work in Ref . |117| . 



Minimal model 

We now set j3 = 0, in which discussed earlier and as checked numerically, the model 

exactly reduces to a collection of one-dimensional FKPP-likc models. Hence, in order to compute 
two-point correlation functions, it is enough to evolve two realizations of the corresponding FIP 
model with the constraint that all bins with k < k&b be identical between the two realizations, 
and the bins k > &A& be completely independent. Alternatively, we could also generate one 
single realization and evolve it for iA& time steps, replicate it at time £a&, and then evolve the 
two replicas completely independently of each other. The difference between these two possible 
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Figure 5.5: The same as in Fig. |5.4|but for (3 = 2 
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Figure 5.6: function of At — t — t&b hi the minimal model with /3 — for N = 10 10 . 

We display the results obtained within the model in which the realizations decorrelate in the bins 
k > fcAh (labelled "Monte Carlo"), and within the model in which the decorrelation is complete 



after time t^ (labelled "without correlations"). The theoretical curves use Eq. (5.401 with the two 
possible choices for the front size L. Inset: The same, as a function of 1/At in order to highlight 
the small-Ai region where, as expected, important differences appear between the models. 
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Figure 5.7: The same as in Fig. 5.6 for N = 10 . All curves coincide almost perfectly. 



implementations of the minimal model cannot be accounted for in our analytical calculations, thus 
the differences that we shall find numerically will give an indication of the model uncertainty. This 
time, our aim is essentially to check our analytical formulas, thus we will pick very large values 
of N, even if they appear to be unphysical in the QCD context since they would correspond to 
exponentially small values of the strong coupling consta nt a s . 

The parameters of the model are obtained from Eq. (5.201 with a — 1, (3 — and dt — 10~ 2 : 



7c 



1.0136- 



v = 2.6817- 



u"(jc) = 2.6098- ■ 



(5.55) 



These values are close to 1, e and e respectively, which would be the correct parameters if dt were 
infinitesimal, in which case u (7) = e> (see Eq. flgggp ). 

The numerical results are shown in Fig. |5.6| for N = 10 10 with the two versions of the model 
(we generated about 10 5 realizations), and compared with the analytical predictions. We test the 
two possible choices for the size L of the front: Either L is a constant, which from our previous 
experience with FKPP traveling waves |114|, we set to 



L = — \nN + — In In AT - — 

7c 7c 7c 



(see e.g. Eq. ( 4.71| ), or it is ^-dependent, namely 

1 



7c 



]nN + S- 1.4. 



(5.56) 



(5.57) 



The numerical constants, which are not determined in our theory, were chosen empirically so that 
they properly describe all numerical data for N > 10 10 . In the first case, Eq. (5.40 1 is used. In the 



second case, Eq. (5.371 is integrated numerically over t and <5. We see that the agreement between 



the numerical calculation and the analytical predictions is good, except maybe for very small 
values of At where the calculations are not expected to be accurate. Indeed, for the same values of 
At, we also see in Fig. |5.6| a sizable discreapancy between the two versions of the minimal model. 
The calculations for N = 10 50 are shown in Fig. 5.7 The numerical results and the theoretical 



expectations (Eq. (5.40 1) coincide almost perfectly. 

For larger and more realistic values of a s , the persistence of the correlations is still seen in 
the numerical simulations, but some parameters should be modified in the analytical expressions 
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Figure 5.8: Comparison of a numerical Monte Carlo simulation and our analytical formula. The 
constant in the parameter L (see the text) which should be equal to ln(l/as)/7o for very small a s , 
has been shifted by a phcnomcnological constant. Once this is done, we get a very good agreement 
between the two calculations. 
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Figure 5.9: Numerical check of the scaling (5.46). The curves for the different values of N are 
very close together for N > 10 10 , but the scaling seems to break down for low values of N (see 
the curve for N = 100), as expected. 
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and tuned to account for our lack of understanding of subleading corrections important for finite 
ln(l/ag). We show such a calculation for a s = 0.1 in Fig. 5.8 compared to a variant of Eq. (5.401. 



Finally, we check that the scaling in Eq. (5.461 is well reproduced by the numerical data. The 
Monte-Carlo simulations are shown in Fig. |5.9[ plotted in the appropriate scaling variables. The 
diffusion constant of a single wave front D as well as the velocity v are measured from the same 
data. We see that all curves nicely superimpose for N > 10 10 (we show data for values of ./V 
as large as 10 80 ), while there are clear deviations for smaller N (see the curve for N = 100), as 
expected. 
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Chapter 6 

Phenomenological applications 



In this chapter, we review the phenomenological consequences of the results obtained from the cor- 
respondence with statistical physics. We derive new properties of the QCD scattering amplitudes 
and discuss their impact on phenomenology. 

Contents 

6.1 Dipole models and geometric scaling 1821 

6.2 Diffusive scaling 



As was stated in the Introduction, the initially unplanned opportunity to collect data in the 
high-energy regime of deep-inelastic scattering at HERA triggered a renewed interest in small- 
x physics among phenomenologists. The major discoveries in this regime is the (unexpected) 
important fraction of diffractive events, and a new scaling, geometric scaling, featured by total 
(and even semi-inclusive) cross-sections (see Fig. |1. 1| in the Introduction). 

In order to deal theoretically with the small- a; regime, one needs new factorization theorems in 
order to single out the elements of the cross-sections that are computable in perturbation theory. 



High-energy, also called fc^-factorization 124 ■ 126 , is the appropriate tool. A practical way to 



implement /c^-factorization is the color dipole model presented in Chap. [2] 



6.1 Dipole models and geometric scaling 

The main observable measured at HERA is the proton structure function F^. It is proportional to 
the sum of the virtual photon-proton cross-section for a transversely and longitudinally polarized 
photon respectively. 

A bare photon has no hadronic interactions, since it does not carry any color charge. However, 
it may easily fluctuate into a quark-antiquark pair, overall color-neutral, thus forming a color 
dipole. Subsequently, these dipoles interact with the target proton. This picture is represented by 
the following equations: 

&t,l{x, Q 2 ) = J dzd 2 r\9 TiL {z,r 1 Q 2 )\ 2 a dipolc (x, r). 

Here, ot,l & re the photon-proton cross-sections for transversly and longitudinally polarized virtual 
photons. ^t,l are light-cone wavefunctions for 7*, computable within QED (see, e.g., Ref. [28] 
for explicit expressions to lowest order in a em ). Furthermore, Odipoio^, r) is the cross-section 
for dipole-proton scattering (for a dipole of transverse size r), and encodes all the information 
about hadronic interactions (including unitarization effects). This cross-section is related to the 
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amplitude A discussed so far by an integration over the impact parameter. (Actually, A was the 
forward elastic amplitude; the optical theorem relates it to the total cross-section). 
In Ref. (28 29 , the dipole cross-section was modeled as 



a d ipoie(x,r) = ^(l-e-^W/*), (6.2) 

where o~o is a hadronic cross-section: It stems from the integration over the impact parameter, 
when the impact parameter dependence is supposed to be uniform over a disk of radius ~ y/o^. 
Q s (x) plays the role of the saturation momentum, parametrized as Q 2 s (x) — (xq/x) x x 1 GeV 2 . 
Note that, by construction, this cross-section only depends on the combined variable r 2 Q 2 s (x) 
instead of r and x separately. This property is transmitted to the measured photon cross-sections 
ct,l(x,Q 2 ), which then depend on Q 2 /Q 2 (x) only (this scaling is slightly violated by the masses 
of the quarks). This is geometric scaling, predicted to be a feature of the solutions to the BK 
equation at large rapidity. 

Historically, geometric scaling was discovered first in the data (see Ref. |30|), after Golec- 
Biernat and Wiisthoff (GBW) had written down their model: The latter happened to feature this 
scaling (up to small violations induced by the quark masses). There was no apparent need for 
finite rapidity scaling violations in the first HERA data. However, later analysis revealed that a 
significant amount of explicit scaling violations in the dipole cross-section, predicted by the BK 
equation, were actually required by more accurate data. 

A now popular model that describes the HERA data in a way that takes a better account of 



the subasymptotics, beyond the GBW model, was formulated in Ref. 127 . The dipole scattering 
cross-section reads (Jdipo\e(x,r) = 2irR 2 Af(y,rQ s ), with 



! fOT rQs " 2 ' (6.3) 

1 _ c -aln 2 (brQ s ) fof r Q s > 2 , 

where Q s = Q s (x) = (x /x) x / 2 GeV. The expression for the cross-section for r small compared 



to 2/Q s corresponds to the solution of the BK equation (compare to Eq. (4.22) with the help of 



Tab. 2.1 1, in which we substituted uj(^ c ) = ^(ic) aid u/'(7 c ) = uj"(^ c ) by the parameters A and 
k that we subsequently fit to the data. The expression in the second line also has the correct 
functional form for r» 2jQ a , as obtained by solving the BK equation [33] . This is strictly valid 
only to leading-order accuracy, but here it is used merely as a convenient interpolation towards 
the 'black disk' limit M = 1. (The details of this interpolation are unimportant for the calculation 
of <r 7 »p.) The coefficients a and b are determined uniquely from the condition that J\f(rQ s , Y) and 



its slope be continuous at rQ s = 2. The overall factor Af in the first line of Eq. (6.3 1 is ambiguous, 
reflecting an ambiguity in the definition of Q s . This model fits well all HERA data for structure 
functions, in the range x < 10~ 2 . All details may be found in Ref. |l27| . 

The model explicitely breaks geometric scaling. However, effectively, geometric scaling remains 
a fairly good symmetry of the model, as required by the data. The small finite-rapidity scaling 
violations are needed to describe accurately the high-precision HERA data. 



The model may also accomodate less inclusive observables, such as diffraction 128 . It has 
been improved recently by including heavy quarks |129| (The crucial need for taking account 
of the charm quark was emphasized in Ref. |130| ). An impact-parameter dependence was also 



introduced 131 ■ 133 that was already missing in the GBW model. 



The range of validity of dipole models has been re-examined recently 1 134 



6.2 Diffusive scaling 

At still higher energies, according to the discussion of Chap. [4] one expects the saturation scale 
to acquire a dispersion from event to event that scales with the rapidity like \fm) when rapidity 
increases. Although this dispersion is not an observable since there is no way to measure the 
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saturation scale of an individual event, it manifests itself in the total cross-section in the form of 
a new scaling, different from geometric scaling. 

The physical amplitude for the scattering of a dipole of size r off some target is given by the 
average of all realizations of the evolution at a given y: 



A(y,r) = (T(r))\ 



v 



(6.4) 



For large enough rapidities and small enough a s , these realizations are exponentially decaying 
fronts in the variable p = ln(l/r 2 ), fully characterized by a stochastic saturation scale, or rather 
its logarithm p s = lnQ 2 (y). For the purpose of the present discussion, it may be approximated in 
the same way as in Eq. (|5.2|, namely 



T(p) = 6( Ps -p) + 6(p - p s )e-^P-<>°\ 



(6.5) 



The statistics of p s is given by Eqs. ( |4. 64 1,(4.66 ) (up to the replacements suggested in Tab. 2.1 
to go from a generic reaction-diffusion to QCD). At ultrahigh energies (and very small a s ), it is 
essentially a Gaussian centered at 



(Ps 



tt 2 7 c u/'(7 c ) 



and of variance 



2(ln(l/a 2 ) + 31nln(l/a 2 )) z 
TrV'frc) 



ay 



(pi) ~ (p.) 



31n 3 (l/a 2 ) 



ay. 



(6.6) 



(6.7) 



The scattering amplitude may be expressed by the simple formula 



i 



dp s T(p)\ y exp 



(Ps - (Ps)f 



ay/2n J ■ " * V 2<j2 
The most remarkable feature of this amplitude is the scaling form for A that it yields: 



(Ps(y)} 



sjay/\n\l/al) t 



(6.9) 



This equation may be obtained by performing the integration in Eq. (6.8) after the replacement 



of T by its approximation (6.5 1. This scaling obviously violates geometric scaling: If the latter 
scaling were satisfied, then A would be a function of p — (p s (y)) only. 

In Ref. 40 , Mueller and Shoshi had already noted that geometric scaling had to be violated 



beyond the BK equation. However, the square root in the denominator of the scaling variable in 



Eq. (6.9 1 was missing because their approach was relying on mean field throughout, thus missing 



the stochastic nature of the evolution. 

This new scaling is a firm prediction of the correspondence with statistical physics. However, 
it may not be tested at particle colliders in a simple way. Let us work out the order of magnitude 
of the rapidity needed for the different effects (saturation, geometric scaling, diffusive scaling) to 
show up. The rapidity that is needed to reach saturation is roughly 



J/BFKL 



ln(l/a 2 ; 
au){\) 



(6.10) 



The BK picture is expected to be valid until the asymptotic exponential shape of the front has 
diffused down to the point where the amplitude becomes of the order of a\. This additional 



rapidity needed to get to the regime of geometric scaling is thus given by Eq. (4.371 once the 
appropriate replacements have been done 



1 



2/BK 



2aw"(7 c ) [ 7. 



ln(l/a 2 



(6.11) 
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and finally, the effect of the fluctuations of the saturation scale gets important at the rapidity 

31n 3 (l/ap , . 

Muct ~ - 3 „, r - (6.12 

The relevant parameters in QCD are deduced from the BFKL kernel. They read 

7 C = 0.627549, lo{ 1c ) = 3.0645, u/'( 7c ) = 48.5176. (6.13) 

For some realistic strong coupling constant, a s ~ 0.2, we get 

J/bfkl ~ 6.07879 , y BK ~ 1.41965 , yfiuct ~ 0.348244. (6.14) 

Given that rapidities in the small- a; regime at HERA were of the order of 10, and will be of the 
order of 15 at the LHC, these figures indicate that we may observe these effects. However, there 
are many criticism to these naive estimates. 

First, the values of the rapidity that delimitate the different regimes are largely underestimated 
given that they rely on the leading-order BFKL kernel, which predicts a much too large growth 
of the cross-section with the rapidity and a too fast diffusion (see the large value of lo"(^ c )). 
Already the effect of the running coupling, which should be taken into account in any detailed 
phenomenological study, is expected to still reduce the effects of the fluctuations 1 135] . 

Next, one also has to keep in mind that the former estimates should only hold for very small 
values of a s , such that \al/a 2 s 3> 1 which is certainly not true in real-life QCD. Note that 
asymptotically, one should have yfl uc t S> dbk 3> J/bfkl- The fact that the order is inverted means 
that the quantitative results obtained within the phenomenological model for front propagation 
should not be trusted for values of a s as "large" as 0.2. 

Nevertheless, the effect of diffusive scaling (i.e. of the event-by-event fluctuations of the sat- 
uration scale) on observables has already been investigated in some detail by several groups. 
Diffractive amplitudes were studied in Ref. [136] . The ratio of the gluon distribution in a nucleus 
to the same quantity in a proton was computed in Ref. |137 . 
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We have reviewed a peculiar way of viewing high-energy scattering in QCD, based on the physics 
of the parton model, and its strong similarities with reaction-diffusion processes (Chap. [2|. The 
correspondence is best summarized in the mapping of Tab. |2.1| We have seen that the equations 
that describe the dynamics of these processes are in the universality class of the stochastic FKPP 
equation, and admit traveling-wave solutions whose features are likely to be universal, in such a 
way that a study of simple reaction-diffusion-like models may lead to exact asymptotic results also 
for QCD scattering amplitudes. Understanding the very mechanism of traveling wave formation 
and front propagation was crucial to see how the universality may come about (see Chap.[4|. 

In zero-dimensional stochastic models, we could perform exact calculations and get analytical 
results within different formulations (Chap. [3j. We understood that analyzing the structure of 
single events was technically much simpler if one wants to get leading orders at large N {= 1/af), 
since in individual realizations, one may factorize the fluctuating part from the nonlinear effects. 
Thanks to this observation, in one-dimensional models which admit realizations in the form of 
stochastic traveling waves, we could also get precise analytical results on the form and shape of the 
traveling waves, which are presumably exact asymptotically (Chap. [4]). Universality enables one to 
make statements on the form of the QCD scattering amplitudes at very high energies. Appropriate 
extensions of the relevant statistical models which incorporate an additional dimension lead to 
predictions for the correlations in the transverse plane (Chap. [5j. Some of these results turn into 
firm phenomenological predictions (Chap. [6]), which however do not seem to be testable at colliders 
in the near future. Nevertheless, getting new analytical results for QCD in some limit is always 
an interesting achievement, given the complexity of the theory. Furthermore, while our analytical 
results only apply for exponentially small a s (ln(l/o!g) 3> 1), the picture itself should be valid in 
the whole perturbative range, namely for Oj<1. 



Prospects. There are still many open questions. On the statistical physics side, the statistics 
of the front position that we have found has not been derived rigorously, but rather guessed, and 
rely on many quite ad hoc conjectures. We got confidence on the validity of our conjectures on 
the basis of numerical simulations. Moreover, although we expect universality up to corrections 
of order 1/N (that is to say 0(a 2 s ) in QCD), we could only get analytical expressions relative to 
the cumulants of the position of the front for the first terms in an expansion in powers of 1/ In N, 
which extremely large values of iV (small a s ) to be valid. But on a more general footing, the 
sFKPP equation seems to describe many physical, chemical or biological problems (in particular 
population evolution with selection in evolutionary biology). We have also found recently an 
explicit analogy with the theory of spin glasses |138 139] . This large universality is maybe the 
strongest incentive to try and find more accurate solutions to that kind of equations. 

On the QCD side, the correspondence with reaction-diffusion processes strongly relies on the 
assumption that there is saturation of some form of the quark and gluon densities in the hadronic 
wave functions. While this is a reasonable guess that few experts would challenge, it is clear that 
we cannot consider that the problem is solved before the saturation mechanism at work in QCD 
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has been exhibited. QCD is formulated as a quantum field theory. To see the similarity with 
reaction-diffusion, we basically needed to translate it into the parton model first. It would be 
better to recover the results of Chap. [4] (and hopefully get more) directly from field theory [59], as 
one could do it in the zero-dimensional model introduced in Chap. [3] This requires to understand 
the strong field regime of field theory. This is an exciting challenge for both particle physicists 
and statistical physicists. 



Let us finally state our personal prospects in the field. First, we wish to go back to the 
simple Balitsky-Kovchegov equation and study in more detail the properties of its solutions. It 
is important for phenomenology, since it seems that in the range of energy that may be reached 
at experiments, effects described by more advanced equations (incorporating genuine saturation 
effects, as discussed at length in previous chapters) are likely to be negligible. Interestingly, since 
the BK equation also represents the statistical properties of the tip of a random walk (see Chap. [4] 
and the recent paper 107| ), some fine properties of scattering amplitudes may be inferred from 
the study of the latter. Work is in progress in this direction in collaboration with Al Mueller. 

On the pure statistical physics side, we wish to pursue the study of simple models like the ones 
presented in Chap. [3j which could be of some interest in the interdisciplinary field of population 
evolution studies. 

Last, if the formalism of the dipole model on which relies most our work in QCD seems well- 
suited for electron-proton or nucleus high-energy scattering, most of the experimental data which 
will become available in the next decade are about proton and nucleus interactions at the Large 
Hadron Collider (LHC). The new challenge to phenomenologists is to formulate and compute 
observables in this context. It seems that quadrupoles play an important role for all interesting 
observables, see e.g. Ref. [140] . Recently, we made a first step in the direction of computing the 



evolution of such objects with the energy 141 and we will pursue in this promising direction. 



An interesting theoretical question in the continuation of our work would be, for example, if the 
correlations computed in Chap. [5] would show up in the evolution of these quadrupoles and hence 
in the corresponding observables measured at the LHC. 
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